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1.0 SUMMARY 


A finite difference method for solving the unsteady transonic flow about 
harmonically oscillating wings is investigated. The procedure is based on separating the 
velocity potential into steady and unsteady parts and linearizing the resulting unsteady 
differential equation for small disturbances. The differential equation for the unsteady 
velocity potential is linear with spatially varying coefficients and with the time 
variable eliminated by assuming harmonic motion. 

The work of this report is a direct extension of earlier studies and includes 
correlation with experimental results, t3q>ical section flutter calculations for a NACA 
64 AO 10 airfoil, and some further theoretical development. 

The main results of this study are as follows: 

1. A study is made of the representation of the shock motion associated with an 
oscillating airfoil. The effects of the shock motion and the resulting pressure pulse 
are shown to be included in the pressure distributions and generalized forces as 
calculated using the harmonic method. 

2. A typical section, two-degree-of-freedom flutter analysis of a NACA 64A010 airfoil 
for Mach numbers from 0.75 to 0.90 is presented. These studies show a sharp 
transonic bucket in one case and abrupt changes in instability modes. 

3. Analytical and experimental pressure distributions for the NACA 64A010 airfoil 
are compared at Mach numbers of 0.75, 0.80 and 0.842. The results are presented 
for plunge and pitch modes and several reduced frequencies. 

4. An alternating direction method is derived for the harmonic method. This 
formulation appears to provide a practical procedure using relaxation solution 
techniques which do not have solution convergence problems. 



2.0 INTRODUCTION 


Aerodynamic flutter may occur in any flight regime, but its presence at transonic 
speeds is of special interest to the flutter analyst; first because the critical flutter speed 
often occurs in this flow regime, and second, because of the difficulty in calculating 
imsteady transonic air forces. The flutter boundary, when plotted in terms of dynamic 
pressure versus Mach number, displays a characteristic "bucket” in the transonic 
regime which may well result in minimum flutter margins. This drop in dynamic 
pressure is partially due to the increase in magnitude of the pressure coefficients as 
predicted by linear theory. An additional reduction in this boundary is due to the 
presence of shocks over the airfoil. The existence of shocks significantly affect the 
aerodynamic forces in terms of magnitude, distribution, and phase angle with respect to 
surface motion, all of which are important in the evaluation of flutter. The degrading of 
the flutter boundary by transonic flow characteristics appears to be particularly 
important for wings with supercritical airfoil sections. This point is discussed by Ashley 
in reference 1. 

The problem in calculating unsteady transonic flow stems from the presence of 
moving shocks and the resulting large velocity gradients. This necessitates the 
consideration of a nonlinear formulation. Also, the flow is "mixed” in the sense that 
both subsonic and supersonic regions exist and thus in adjacent regions in the flow field 
the governing differential equations may be of distinctly different character. 

The current means for obtaining reliable flutter characteristics in the transonic 
region is through wind tunnel experiments. This expensive and time consuming process 
makes it highly desirable to have efficient analytic procedures available for the 
evaluation of transonic flutter so that the complete flutter characteristics can be studied 
during all phases of vehicle design and development. 

In recent years, the steady-state transonic flow about a lifting surface of finite 
thickness has been satisfactorily handled using a finite difference formulation of the 
nonlinear differential equation and relaxation techniques to solve the resulting set of 
simultaneous equations. The results have correlated well with experimental data. The 
procedure is time consuming but appears practical and compatible with generally 
available resources. 

The unsteady problem is not, as yet, capable of being handled in as practical a 
fashion. First of all, it is a problem of considerably more scope, as air forces must be 
generated for mode shapes and reduced frequencies as well as for the Mach number of 
the steady problem. Thus, the number of cases to be computed is significantly larger. 
Secondly, since the additional coordinate of time is included, it has meant a 
significantly more complicated formulation - either imposing a time integration on top 
of the spatial relaxation solution, or assuming harmonic motion and working with the 
resulting complex spatial differential equations. 
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It is possible to calculate the pressure distributions for unsteady motion using the 
complete nonviscous, compressible, time dependent differential equations. This has been 
done, for example, by Magnus and Yoshihara and by Beam and Warming (references 2 
and 3). The results are obtained using a finite difference representation of time as well 
as space derivatives. However, the cost of doing these analyses make this procedure 
impractical for flutter analyses. This is also a nonlinear analysis which increases the 
number of cases to be considered, since superposition of solutions is not possible. 

A reasonable simplification is achieved by assuming transonic small perturbation 
from a uniform parallel flow. The resulting time dependent equations have recently 
been used by Ballhaus and Goorjian (ref. 4) to generate a program called LTRAN2 for 
analyzing the two-dimensional, low-frequency problem. LTRAN2, which uses an implicit 
finite difference algorithm to integrate the nonlinear differential equation, is very 
efficient. This program has been extended by Rizzetta, Yoshihara and Borland to an 
EXTRAN2 version which includes high frequency terms in the differential equations 
and boundary conditions (ref. 5), viscous interactions (ref. 6 and 7), and three 
dimensions (ref. 5). 

An alternative to this approach is achieved by separating the velocity potential 
into steady and harmonically varying unsteady parts. The result is the usual nonlinear 
differential equation for the steady velocity potential, and a linear equation for the 
unsteady velocity potential, with spatially varying coefficients which are functions of 
the steady velocity potential. This is the procedure to which this work is directed. The 
basic derivation is provided by Ehlers in reference 8, with the results of additional 
development and exploration described in references 9 through 13. The use of the steady 
state potential in the unsteady differential equations introduces the steady state flow 
properties into the propagation of the disturbances from the oscillations of the wing 
including the shocks. It has been sometimes considered that the solution of the complex 
differential equation for harmonic motion did not include the effect of shock motion on 
the unsteady forces. It will be shown that the frequency domain solution does indeed 
implicitly take into account the lift forces and moments produced by the moving shock 
wave. 

In the following, calculations of the pressure distributions for several Mach 
numbers and reduced frequencies will be presented and the results compared with the 
corresponding experimental measurements made on an NACA 64A010 airfoil at the 
NASA Ames Research Center by Davis and Malcolm (ref. 14). The frequency domain 
method will also be applied to a two-dimensional flutter analysis which will illustrate 
the transonic effects of the moving shock on the characteristics of flutter. These flutter 
results are directly comparable to the results of a linear analysis by Isogai (ref. 15) 
using air forces which are the first harmonic components from a time integration of the 
transonic small perturbation equation. 

A two step method, which uses the linearized solution for a flat plate to reduce 
significantly the cost of computing unsteady air loads for a thick airfoil by the 
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frequency domain technique, will be described. An ADI method of solving the complex 
differential equation for the potential will be presented which does not have the 
frequency limitation of the conventional relaxation method. This method may reduce 
the amount of computer resources for the three-dimensional problem, making three- 
dimensional transonic analysis practical. 

The authors would like to acknowledge the valuable assistance of Dr. Donald 
Rizzetta for calculating all the EXTRAN2 results. 



3.0 ABBREVIATIONS AND SYMBOLS 


a Streamwise dimension of mesh region; also amplitude of pitching or plunging 

modes, also elastic axis position (see fig. 42) 

an Fourier coefficients for sin(n0l) 

A Generalized force matrix 

b Root semichord of wing or semichord of airfoil; also vertical dimension of 

mesh region 

bn Fourier coefficients for cos(n0i) 

Cp Pressure coefficient, (p - po)/(l/2 poUo^) where p is the local pressure, po the 

freestream static pressure, and pQ the freestream air density 

Cl Lift coefficient based on chord 

Cl« Lift coefficient based on chord for pitching motion 

CLh Lift coefficient based on chord for plunging motion 

Cm Moment coefficient based on chord 

Cm« Moment coefficient based on chord for pitching motion 

CMh Moment coefficient based on chord for plunging motion 

f(x,y,t) Instantaneous wing shape defined by zq = Sf(x,y,t); also shock location 

f Frequency in Hertz 

fO Undisturbed wing or airfoil shape 

fl Unsteady contribution to wing or airfoil shape 

g Structural damping 

h = hei<^f, airfoil plunging motion, positive down 

I x,y subscripts and indices for points in the mesh 

i VT 
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Pitching moment of inertia 


la 

k Reduced frequency based on semichord, 27rfb/U. Same as (o, 

K Transonic parameter, (l-M2)/(M2e); also generalized stiffness matrix in 

section 7 

le Leading edge 

L Total lift; positive upward 

M Freestream Mach number; also generalized mass matrix 

My Total moment about elastic axis 

n (nx, ny, nt), unit normal vector to shock 

q o)2/€ - i(t){y - D^n 

^xx 

T(x Dimensionless radius of gyration about center of gravity 

Sa Static unbalance about elastic axis, positive aft of elastic axis (see fig. 42) 

t Time in units of b/U 

te Trailing edge 

u K - (y + 1) <^0x 

U, V Freestream velocity 

xQj yO Physical coordinates, made dimensionless with the root semichord 
X, y Scaled coordinates (xq, />tyQ) for the two-dimensional problem 
Xq Steady chordwise shock location 

Xi Complex amplitude of shock oscillation 

Xm Magnitude of Xi 

Xa Static unbalance about elastic axis, positive for center of gravity aft of elastic 

axis 

a Angle of attack: also pitching motion about elastic axis 
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I 


am Mean angle of attack 

“ [^] section 5 and Appendix A; e-i^^t in section 8.2 and Appendix D 
P2 = <ux> 

y Ratio of specific heats for air 

ACp Jump in pressure coefficient across airfoil or wake 

A<p^ Jump in at plane of wing or vortex wake 

A(p^te Jump in <pi at wing trailing edge 

6 Thickness ratio, also finite difference operator 

€ (8/M)2/3 

\1 a>M/(l - m2) 


fjL Scale factor on yo, = Sl/3 m2/3j also airfoil to air mass ratio, m/(7rpb2), 

where m is the sectional mass per unit span and p is the freestream density 

V Phase of shock amplitude X\ 

0i = (ot V 

p Free stream density 

(p Complete, scaled perturbation velocity potential; also used for the unsteady 

potential in finite difference equations with subscripts 

Steady scaled perturbation velocity potential 

Unsteady scaled perturbation velocity potential 

V?2 = 

Flat plate complex potential; also linear time dependent potential 
\{f Time dependent potential 

0 ) Reduced frequency in radians; same as k 
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oia 

wh 

[] 

< > 

A 

[] 

{} 

[ 1 ] 

F 


Pitch mode resonant frequency, radians per second 
Plunge mode resonant frequency, radians per second 
Notation 

Denotes jump in quantity across shock 
Denotes mean value of quantity at shock 

Denotes jump in quantity across airfoil; also increment of lift and moment 
coefficient from shock motion in section 5.5 

Matrix Notation 

Rectangular matrix 

Column matrix 

Unit matrix 

Subscripts 

Flutter condition 


1 1 II III I III II II I 



4.0 FORMULATION AND SOLUTION 


Since the mathematical derivation of the method for the solution of the unsteady 
velocity potential for the flow about a harmonically oscillating wing is presented in 
reference 8, the discussion here will be limited to a brief outline of the procedure for 
two dimensions. The complete nonlinear differential equation was simplified by 
assuming the flow to be a small perturbation from a uniform stream near the speed of 
sound. The resulting equation for unsteady flow is 

[K - (7 - 1 V, - (7 + 1 )». J + Vyy - (2»xt + %) / " = 0 M l) 

where K = (1 - M2) / (M2 e), M is the freestream Mach number of velocity Uq in the 
x-direction, x and y are made dimensionless to the semichord b of the airfoil and the 
time t to the ratio b/UQ. With the airfoil shape as a function of time defined by the 
relation 

Yq = 5f(x,t) 

the linearized boundary condition becomes 

(Py = fx(x,t) + f^(x,t) (4.2) 

The quantity 6 is associated with properties of the airfoil (such as maximum 
thickness ratio, camber, or maximum angle of attack) and is assumed to be small. The 
coordinate y is scaled to the dimensionless physical coordinate yq according to 

y = 5^/^ 

and € is given in terms of 8 by 
e = (5/M)^/^ 

The pressure coefficient is found from the relation 

The preceding differential equation is simplified by assuming harmonic motion 
and by assuming the velocity potential to be separable into a steady-state potential and 
a potential representing the unsteady effects. We write for a perturbation velocity 
potential 

= (/>Q(x,y) + iP|(x,y)e^^^ (4.3) 

and for the body shape 

yo = Sf(x,t) = 6 [fo(x) + fi(x)e^‘^t]. 

Since the steady-state terms must satisfy the boundary conditions and the 
differential equation in the absence of oscillations, we obtain 
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with 


v?Oy = fo(x), y = 0, -1 < X < 1 

On the assumption that the oscillations are small and products of <p\ may be 
neglected, equations (4.1) and (4.2) with the aid of equations (4.4) and (4.5) yield 


(4.5) 


where 




+ - (2iw/e)<^2^ + = 0 


q = oP"! € - io;(7 - 1 

subject to the wing boundary conditions 


(4.6) 


(^2 = fj + icofj(x), y = 0, -1 < X < 1 


(4.7) 


A computer program for solving the steady state transonic flow about lifting 
airfoils based on equations (4.4) and (4.5) was developed by Cole, Murman, and Krupp 
(refs. 16 and 17). For some calculations of unsteady flow, the latest version called 
TSFOIL (ref. 18) was used to compute the coefficients for the difference equations. 
Steady-state solutions from the time dependent method of Rizzetta and Chin (ref. 19), 
called EXTRAN2, to solve equation (4.1), were also used. An additional improvement on 
locating the shock is also available both in TSFOIL and EXTRAN2. For airfoils at high 
Mach numbers and angles of attack, TSFOIL and EXTRAN2 predict shock positions 
considerably aft of experimental measurements. To overcome this difficulty, Jou and 
Murman (ref. 20) developed a phenomenological model for the displacement thickness 
effects of shock wave boundary layer interactions. A wedge (or ramp) is introduced 
behind the shock to simulate the thickening of the boundary layer. Unsteady flow 
examples are presented with the steady state potential computed both with and without 
the viscous ramp. 


The similarity of the unsteady differential equation to the steady-state equation 
suggests that the method of Cole, Murman, and Krupp should be an effective way to 
solve equation (4.6) for the unsteady potential <p\. Note that equation (4.6) is of mixed 
type, being elliptic or h 3 q)erbolic whenever equation (4.4) is elliptic or h 3 q>erbolic. 
Central differencing was used at all points for the y derivative and all subsonic or 
elliptic points for the x derivatives. Backward (or upstream) differences were used for 
the X derivatives at all hyperbolic points. 

The boundary condition that the pressure be continuous across the wake from the 
trailing edge was found in terms of the jump in potential to be 

e (4.8) 

where is the jump in the potential at x = xte just downstream of the trailing edge 

and is determined to satisfy the Kutta condition that the jump in pressure vanish at the 
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trailing edge. The quantity A^i is also used in the difference formulation for the 
derivative ^lyy to satisfy continuity of normal flow across the trailing edge wake. 

Fot the set of difference equations to be determinate, the boundary conditions on 
the outer edges of the mesh must be specified. In the original unsteady formulation, 
these boundary conditions were derived from asymptotic integral relations in a manner 
parallel to that used by Klunker (ref. 21) for steady flow. A later formulation (ref. 10) 
applies an outgoing plane wave boundary condition to the outer edges of the mesh. This 
boundary condition is numerically simpler to apply and is equivalent to the 1st order 
nonreflecting boundary conditions derived by Engquist and Majda (ref. 22). 

The preferred numerical approach to solving the resulting large-order set of 
difference equations is a relaxation procedure, which permits the calculation to be made 
as a sequence of relatively small problems. However, as discussed in proceeding NASA 
reports by the authors (refs, 9 and 10), a significant problem of convergence with the 
relaxation procedure was encountered that severely limits the range of Mach number 
and reduced frequency for which solutions may be obtained. Accordingly, an out-of-core 
solver (ref. 23) was developed to solve the complete set of difference equations 
simultaneously, which for two-dimensional flow is fairly efficient. Recently, an ADI 
method based on the time dependent method to obtain frequency domain solutions was 
derived and preliminary results indicate no frequency limitation on convergence. This 
method is derived in Appendix D and the results discussed in section 8.2. 



5.0 INVESTIGATION OF SHOCK WAVE MOTION ON AN 
OSCILLATING AIRFOIL 


5.1 DESCRIPTION OF BASIC SHOCK MOTIONS AND EFFECTS ON FLUTTER 

Tijdeman (ref. 24) has experimentally observed and measured the motion of shocks 
on an airfoil with an oscillating flap. He has found three t 3 rpes of shock motions. Type A 
is sinusoidal shock-wave motion. This is characterized by a shock whose position varies 
sinusoidally along the chord as the flap is oscillated sinusoidally. Even for small 
frequencies, there is a significant phase shift between shock positions and flap motion. 

A change in shock strength is also observed with the shock strength attaining its 
maximum at its farthest upstream position. T 3 rpe B is interrupted shock wave motion 
and is characterized by the disappearance of the shock wave during its backward 
motion. Type C is the upstream propagated shock wave. When the flow is just 
supercritical, a shock wave forms on the upper surface as the flap just passes the 
maximum downward deflection. The shock wave moves upstream at first gaining in 
strength, continues to move upstream while decreasing in strength, then leaves the 
airfoil at the leading edge. This phenomenon is repeated periodically. 

Ashley (ref. 1) and McGrew et al (ref. 25) indicate that the predominant factor in 
the failure of linearized theory to predict accurately the forces in transonic flow is the 
presence of shock waves on the upper surface and sometimes the lower surface of the 
wing. The motion of these shocks along the chord can introduce large unsteady 
moments. Ashley (ref. 1) in a simple approximate analysis has shown that this shock 
motion usually destabilizes single degree pitching motion and can have a profound 
effect on flexure-torsion motion. 

It has been frequently stated that the existing frequency domain methods (i.e. the 
methods of Ehlers, in reference 8, and Traci et al., in references 26 and 27) assume a 
fixed shock location, and hence do not yield forces which include the effect of shock 
motions. It was thought that special treatment was required to incorporate these efforts. 
However, even the earliest computations in reference 8, show the complex pressure 
pulse at the steady state shock location which was observed by Tijdeman in his 
experiments (ref. 24). We can find out whether a frequency domain solution represents a 
moving shock by combining the complex pressure Cpi from the frequency domain 
solution with the steady state pressure coefficient CpQ in the form 

By varying the time t we obtain the pressure coefficients shown in figures 1 and 2 
(these results do not contain the shock point operators discussed in section 5) for the 
upper and lower airfoil surfaces. The amplitude of the pitching NACA 64A010 airfoil is 
+2 degrees about a mean angle of attack of -0.21 degrees. The Mach number is 0.8 with 
a reduced frequency of 0.247. From these graphs we see that the shock position does 
indeed move although the amplitude of motion may be much smaller than the 
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amplitude experimentally observed. Since the frequency domain method is linear, only 
type A shock motion is possible from the theory. In sections 5.2, 5.3, and 5.4, we shall 
show that the effect of shock motion is indeed incorporated in the frequency domain 
method used here. It should also be noted the shock motion effects are included in the 
procedure used by Traci et al., in references 26 and 27. 

5.2 EFFECT OF GRID REFINEMENT ABOUT THE STEADY STATE SHOCK ON 

THE SOLUTION 

For strong shocks, the complex pressure pulse did not appear to be accurately 
I' described by the initial grid which was relatively coarse in the neighborhood of the 

shock. Accordingly, the grid was refined in the vicinity of the the shock in an attempt to 
improve the definition of the pressure pulse. Figures 3 and 4 show the jump in pressure 
coefficient across the airfoil for three successively refined grids for M = 0.8, a mean 
angle of attack of one degree, and a reduced frequency of 0.3. This example was chosen 
for study because of the presence of a strong shock which occcurs on only one side of the 
airfoil. Since grid refinement does not define the pulse more accurately it appears that 
the pressure is singular at the steady state shock location. Since this singularity is of 
lower order for the perturbation potential than for the pressure coefficient which 
involves the first derivative with respect to x, we examined the jump in <pi across the 
airfoil for the same three potential distributions. These results are plotted in figures 5 
and 6. It is seen that, as the grid is refined about the steady state shock location, a 
jump in the unsteady potential appears across the shock. 

In Appendix A, the jump conditions across the shock are derived. The analysis 
follows the work of Murman, Hafez, and Riszk (ref. 28), who expanded the shock 
conditions about the steady state location. The condition that the velocity potential be 
continuous across the moving shock is applied at the fixed location of the steady state 
shock, leading to the continuity of (po, the steady state potential. Also, for the unsteady 
potential, the complex amplitude of shock motion is defined by the jump in the unsteady 
potential and the pressure jump across the steady shock.This is given by the relation 

= 0 (5.2.1) 

where Xi is the complex amplitude of the shock motion and denotes the jump in the 
quantity across the steady shock. Hence it appears that the solution with the refined 
grid does indeed have the essential information for calculating the effects of the moving 
shock, since once we find the jump in the amplitude of shock motion is found from 
the foregoing formula. 

The effect of refining the grid on the perturbation potential was further 
investigated by examining three dimensional representations of the potential 
distributions, portions of which are shown in figures 7 to 9. The developing of the 
jump in (PI at the steady state shock location is apparent. Figure 9 exhibits a sharp 
spike away from the airfoil at the shock. On a physical basis, the spike would be 
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expected to occur at the airfoil surface. This is undoubtedly due to some of the points 
along the column of the shock being treated inconsistently as elliptic, while others are 
hyperbolic because of errors introduced from interpolation of the steady state potential 
in refining the grid about the shock. As we shall see, this problem is remedied by the 
introduction of a shock point operator which is described in the next section. 

5.3 EFFECT OF SHOCK POINT OPERATOR ON THE SOLUTION 

To capture steady state shocks accurately, Krupp and Murman (ref. 17) introduced 
a shock point operator to the difference formulation of the transonic small perturbation 
differential equation. A similar shock point operator was derived in reference 10 for 
unsteady harmonic transonic flow and applied to the three dimensional problem. The 
differential equation is written in conservation form and the divergence theorem is 
applied to the derivative terms for an infinitesimal rectangle enclosing points on 
opposite sides of the shock. The derivation is repeated in Appendix B with some 
improvements. Three forms of the difference equation for are presented. The 
original form yields a weighted average of upwind and central differences for at the 
point ij. The second uses a weighted average of central differences for (p\^ at the points 
i and i-1. The third is essentially <p\^ defined by the difference between values of <p\ at 
the right and left sides of the control volume for the shock point operator (see figure B-1 
and the discussion in appendix B). Forms 2 and 3 are identical for equally spaced 
points. 

When the shock point operator is introduced, the solution exhibits a jump in <p\ at 
the shock for the coarse grid as well as for the fine grid. Figures 10 and 11 are three- 
dimensional graphs of the real and imaginary parts of on the xy grid plane for 
M = 0.8, a mean angle of attack of 1®, and the NASA 64A010 airfoil pitching about the 
quarter-chord point at a reduced frequency of k = 0.3. The jump in unsteady potential, 
across the airfoil is given in figure 12 for the same solution using the first form of 
the shock point operator. The jump in <p\ across the shock is most easily seen in the 
imaginary component of in figures 11 and 12. Williams (ref. 29) has shown that the 
perturbation potential at the shock location has a singularity on the downstream side of 
the shock of the form x log x at x = 0. The finite difference solution in figures 11 and 12 
would seem to indicate such a singularity. When the third form of the shock point 
operator is used, however, the jump in <p\ is more pronounced in the real part but there 
appears to be no evidence of the x log x singularity in the real part although the 
imaginary part is essentially the same (see figs. 13 and 14). 

5.4 EFFECT OF SHOCK POINT BOUNDARY CONDITIONS ON THE SOLUTION 

We have seen that for the shock point operator in the presence of strong shocks 
that a jump in <p\ occurs at the steady state shock location. It is also possible to replace 
the shock point operator with the shock boundary conditions in equation (A-23). These 
conditions were coded and applied to the rows near the airfoil where the shock is 
strongest. The shock point operator was applied to the remaining rows. For the 
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boundary conditions to work successfully it was found necessary to use a fitted shock 
instead of the captured shock where the shock effects are smeared over 3 or 4 grid 
spaces. Since we had no steady state program with a shock fitting routine (see, e.g. 
Seebass'et al, ref. 30), we manually constructed fitted shocks for the first 5 rows away 
from the airfoil. Figure 15 shows some examples by the plot of -u where 
u = K- (y+l) ^Ox fhe rows. To obtain a reference result for comparison, we also 

used the same u distribution containing the fitted shock in a conventional solution 
using the shock point operator. The extrapolated jump in <pi across the airfoil is 
compared in figure 16 for the two solutions. The Mach number is 0.8 and the NACA 
64A010 airfoil was pitched about the quarter chord point at a mean angle of attack of 1® 
with a reduced frequency of 0.3. The amplitude of shock motion for the two solutions is 
about the same but the phase differs considerably. The value of Xi computed by 
equation (5.2.1) is 0.769-3.493i for the shock boundary conditions and 1.735-2.090i for 
the shock point operator. For a pitch amplitude of 1/2*^, using the shock boundary 
conditions gives a somewhat larger shock motion of 3.1% of the chord compared with 
2.4% of the chord for the shock point operator. The shock point operator with the 
captured shock gives Xi = 1.62-2.31i with a motion of 2.4% of the chord. For the shock 
point operator, a fitted shock has little effect on the result. 

Since the shock point operator and shock boundary conditions give somewhat 
different results, it seems appropriate to see what the difference is between the jump 
conditions at the shock resulting from the two relations. From equation (A-21) the 
boundary condition that was applied is 

<2iw/e + [<^i] - 0 (5.4.1) 

where [ ] denotes the jump in the quantity across the steady state shock, while < > 
denotes the average value at the shock. Applying the divergence theorem about a 
section of the shock to the differential equation for <pi yields 

[u>Pi^-2We]nx+[‘(’y]ny = 0 

Since ny = 0 and [u ] = 0 for a normal shock, using the identity 

[Wl J = <“> [v-lj ^ <^1,> 1»1 


yields 


■ <2ico/e> [<^l] =0 

for the jump across the shock. The term Ux is seen to be lacking. This is the 
contribution from expanding u on the moving shock about the steady state shock 
position. It appears that at the shock point we can satisfy the difference equation at the 
shock but the shock condition in equation (5.4.1) is not satisfied. If we satisfy the 
correct shock relations of equation (5.4.1) then we cannot satisfy the difference 
equation. 
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The shock point operator is easier to apply since it does not require special fitting 
of the shock and is the form used for all subsequent calculations discussed in this 
document. 

When the shock is weak, there is no distinct jump in <p\ at the steady state shock. 
Figure 17 shows the real part of the jump in across the airfoil for the NACA 64A010 
airfoil at M = 0.8 and zero angle of attack, using the shock point operator and a 
uniform grid about the shock. Even refinement of the grid as in Figure 18 does not 
produce a definite jump in (p\ at the shock. Note that in the vicinity of the shock the 
slope changes rapidly. This change in, the slope produces the* large pressure pulse shown 
in the examples in section 6.0. 

5.5 INCORPORATING THE EFFECT OF THE MOVING SHOCK WAVE IN THE 
CALCULATION OF LIFT AND MOMENT 

Once the solution <pi is found then the amplitude of shock motion is found from 
equation (A- 15); that is 

where f^i = [u] = [K - (7 + D^Oxl- Since the jump in (pi is complex then we write 

where Xm is the magnitude of Xi and v is the phase angle. The shock motion about the 
steady state position is given by 




where a is the amplitude of oscillation of the airfoil motion, here assumed small. Let 
6 1 = B + V then the shock motion is described by 

i0i 


aXj„e 


1 


Figure 19 illustrates the pressure pulse produced on the airfoil by the moving 
shock wave from Tijdeman (ref. 24). The unsteady pressure pulse is seen to be 

[<^0x] H(x-xo-aXj„sin0i) 

where and H(t) is the unit step function which equals 0 for t < 0 

and equals unity for t > 0. Here + and - denote downstream and upstream values. We 
now compute the Fourier coefficients for the pressure pulse. Expressing the Fourier 
series in the form 

00 00 

X) aj^ sin n0 1 + 5!^ cos nd j 
n=l n=0 
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we find the values of an and bn are given by 

an=7 J (-2epoj)H(x-xo-aXj„sin0i)sinn0id0i 

\ = ^J (-2e pO ■ ""O - sin 0 1 ) cos n0 j 

-7T 


We assume x > xq but x - xq < aXm- Then 



The terms of the Fourier expansion of the pressure pulse for the first harmonic become 


aj sin 0 j = 


7T 


- xqX- 

fl-l-r: ) sin(t> + wt) 


.aX. 
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4e 


= Imag ■ 




fX - Xr 


1 _(_J\ 


laX 


m 


Choosing the pressure pulse as defined by 
Imagj(c + id)ei*^^J 



we see that the complex pressure pulse is given by 


c + id 


/ /■‘-■‘o’ 

V \^m , 


The contribution to the lift and moment from this pressure pulse can be integrated 
in closed form. Thus for the lift we have 



where we have used equation (A-15). The 2 factor is needed to make the lift coefficient 
based on chord rather than semichord. Note that the contribution to the lift is 
equivalent to a pressure distribution 

-2e[^i]a5(x-Xo) 
where 8(t) is the Dirac delta function. 


For the moment we have 


iX 


m 


Let t = (x - xo)/aXm> then 


4AC^ = 


7T 


/^XQ+aXjn 

J ^ 

•'xQ-aXjn 
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fx “ Xr 
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dx 




[xq “ Vl - dt 


4ACm= (xo'Xa)ACL 
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To find an estimate of higher harmonics on the lift and moment, we consider the 
second ha^^monic terms, these are 

a = 0 


2e 


b 2 cos 20 j = - 


N 


sin 2 ( sin 


7T 


m ^ 

- 


jcos 2(v + cot) 


= Real 




'x-Xo\ / /x-Xo^ 
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ttX. 
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,aX 
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The complex lift from this term is 
2e 


2ACl = - 


PojXi^ 


7tX 
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2e 
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The moment is seen to he 


2e 






ttX. 
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X - Xr 
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VTT? 
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dt = 0 


aX 


dx 
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^xq - + aXj^tJ t Vl - dt 


= -^[%]XiV/4 = e[^l]a2xi/4 


The lift from the second harmonic is seen to be zero and the moment is second 
order in the amplitude of the airfoil motion compared to first order for the first 
harmonic. Thus the higher harmonics in the shock motion can be neglected. 


We have seen that in the solution a jump in tpi occurs across fairly strong shocks. 
To compute the total lift coefficient we must integrate 

^1 

2Cl = - 2ea J | + icoAiP| + |^A<^jj5 (x - Xq) | dx 


where we have added the contribution from the moving shock wave in the form of a 
delta function distribution. The quantity can be integrated numerically in spite of 
the jump in Aipj. 



Consider the combination 

A 

For X greater than xq, we account for the discontinuity in <p\ in integrating <p\^. Thus 

/ Av?j .dx = -A<^j(-l) + A<^i-A</J2 ''' + A>^2(1) + [A(/jjJ 
= ^'^lte 

Thus the integration of ^<Plx pl'is the contribution from the shock motion gives the 
value of A^i at the trailing edge found from the integration. The lift is then given by 

.1 


2Cl = - 2ea 


^'^Ite' 


/■ 


icoA<Pjdx 


This is the same form obtained for continuous A(pi. A similar analysis for the 
moment coefficient yields 

A<Pite+ / (iojx - l)A>Pidx 
-1 


4Cj^ = - 2ea 


-^a^L 


These two formulas were also obtained by Williams (ref. 29). 


When the trapezoidal rule of integration is applied, no special consideration of the 
shock is necessary. When a Simpson’s rule type of integration is employed, the 
immediate step in (pi at the shock is treated by the trapezoidal rule with the remainder 
of the wing integrated by the Simpson’s rule formula. For the grid distribution 
currently employed, integration by the trapezoidal rule agrees to four places with the 
Simpson’s rule result. This is the method used in all the calculations presented. 
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6.0 CORRELATION OF THEORETICAL CALCULATIONS WITH 
THE EXPERIMENTS OF NASA AMES RESEARCH CENTER 


Calculations by the frequency domain method were made at Mach numbers and 
reduced frequencies which can be compared with the experimental measurements of 
Davis and Malcolm (ref. 14) on a NACA 64A010 airfoil at the NASA Ames Research 
Center. Only the zero mean angle of attack cases were studied since the 4^ results 
appear to contain an excessive amount of separation which the theory is not able to 
take into account. 

6.1 CORRELATION OF STEADY STATE PRESSURE DISTRIBUTIONS 

It is essential that the theoretical steady state pressure distributions also agree 
with the measurements, since the steady state potential is used in the coefficients of the 
equations for the calculation of the unsteady pressures. The steady state potential 
distributions were computed by the program EXTRAN2 described by Rizzetta and Chin 
(ref. 19), although some calculations were made using TSFOIL (ref. 18). For identical 
boundary conditions the two programs yield almost identical results. For a few cases 
imsteady pressures were also computed by EXTRAN2. The first harmonic of the Fourier 
expansion was calculated from the time dependent solution to compare with the complex 
pressures of the frequency domain method and with the corresponding NASA 
experimental data. 

Figure 20 compares the EXTRAN2 steady state pressure distributions on the 
airfoil with the NAS A- Ames measurements at a Mach number of 0.75 and a zero angle 
of attack. Agreement is seen to be quite good. In this example the flow is completely 
subsonic. The NACA 64A010 airfoil used by NASA is not quite symmetric. For the 
calculations included here, the measured airfoil coordinates listed in Davis and Malcolm 
(ref. 14) were used for the boundary conditions, with the leading edge slopes modified by 
Riegel's rule in which the slopes s are replaced by sV 1 + s2 and with some further 
changes to reduce the excessive pressure rise at the leading edge which occurred when 
all the slopes were determined by the coordinates. Riegel’s rule eliminates the 
singularity in the slope at the leading edge of a blunt airfoil. 

Figure 21 compares the steady state pressures computed by three methods for a 
Mach number of 0.8 and zero angle of attack. The dotted curve is the EXTRAN2 results 
with the Riegel’s rule modification. The other two curves are obtained by using TSFOIL 
with the NASA Ames coordinates and two different scalings but without the application 
of RiegeFs rule. The classical Spreiter scaling is that presented in section 4 and follows 
naturally from the procedure of mathematically simplifying the differential equation. 
The Krupp scaling is empirical and leads to the following definitions of K and e: ' 

K= (l - M^) / (m52/3) , 6 = 8'^^^ I 

Although all the unsteady calculations were made with the Spreiter scaling, the 
program is capable of using either parameterization. 


The shock is farthest aft and is stronger for the empirical Krupp scaling than for 
the conventional Spreiter scaling which EXTRAN2 also uses. All subsequent unsteady 
calculations were made with Spreiter scaling. The application of Riegel’s rule improves 
the pressure distribution in the leading edge region of the airfoil but the shock wave is 
weaker than the measured shock. 

Figure 22 shows a graph of the upper and lower pressure distributions for 
M = 0.842 and zero angle of attack. Here the agreement between experiment and theory 
is better, with the shock strength more accurately defined by the calculations than for 
the 0.8 Mach number case. 

6.2 CORRELATION OF UNSTEADY PRESSURE DISTRIBUTIONS 

Using the steady state potential distribution for M = 0.75 represented by the 
pressure distribution in figure 20 to evaluate the coefficients of the unsteady difference 
equations, the complex unsteady pressure distribution was calculated for the NACA 
64A010 airfoil oscillating in pitch about the quarter chord point at a reduced frequency 
of k = 0.20. Figure 23 shows plots of the real and imaginary parts of the jump in 
pressure coefficients across the airfoil for both theory and experiment. Also shown is the 
flat plate, solution for M = 0.75 and k = 0.2. Since the flow is completely subsonic, the 
linearized and the frequency domain solutions are both smoothly varying without a 
large pressure pulse. Both the experimental results and the airfoil calculations show a 
greater pressure amplitude near the leading edge than the linearized or flat plate 
solution. Even for strictly a subsonic case, the influence of the steady flow potential in 
the coefficients of the differential equation has a significant effect on the pressure 
distribution. 

Figures 24 through 31 show the jump in pressure coefficients across the airfoil 
oscillating in pitch about the quarter chord point for M = 0.8, using the steady state 
potential from EXTRAN2 represented by the dotted line pressure in figure 21. In 
figures 24 and 25, the first harmonic from an EXTRAN2 time dependent solution is also 
included with the frequency domain solution for k = 0.101. In the nonlinear EXTRAN2 
solution, the shock pressure pulse is aft of that of the frequency domain solution and is 
closer to the experimental location. In the frequency domain solution, the shock wave 
pressure pulse is located at the steady state shock position. The mean shock location of 
the time dependent solution is generally downstream of the steady state shock position. 

An examination of fibres 24 through 31 indicates that the current linear theory 
locates the shock pressure pulse a distance between 5 and 10% of the chord upstream of 
the measured position. In figures 24 and 25 for k = 0.101, the amplitude of the shock 
pulse from the frequency domain solution is closer to the experiment than the 
EXTRAN2 calculations. For k = 0.202, figures 26 and 27, the magnitude of the pressure 
coefficients from the frequency domain solution compares favorably with experiment, 
particularly in the region between the leading edge and the shock. For k = 0.247, 
figures 28 and 29 compare the pressure jump from both the time dependent method and 
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frequency domain method with experiment. In the leading edge region in figure 28, 
EXTRAN2, which is dependent on amplitude and uses the amplitude of the 
experiment, overpredicts the real component of the pressure, while the frequency 
domain method, which is independent of amplitude, is in better agreement with 
experiment. Aft of the shock, the two theories are fairly close but the real part has a 
sign opposite from the experiment, although the magnitude is small. 

Figures 30 and 31 compare the two theoretical calculations with experiment for 
k = 0.303. Except for the sharp pulse, the two theories compare well for the real 
component of the pressure jump. The frequency domain solution gives a higher peak in 
better agreement with the measurements. The experimental results exhibit greater 
oscillation with position along the chord. The large peak in the real part just upstream 
of the shock pressure pulse at about one-third chord is not predicted by either theory. It 
may be entirely eliminated by removing a single point from the measured data and 
hence may possibly be attributed to scatter in the experimental data. 

The frequency domain solutions for the NACA 64 AO 10 airfoil oscillating in pitch 
about the quarter chord point at a Mach number of 0.842, using the steady state 
potential from figure 22, is given in figures 32 and 33 for a reduced frequency of 0.202. 
For this example, the jump in pressure in the region from the leading edge to the shock 
is fairly accurately predicted by theory. However, the theory yields a large pressure 
pulse at the shock which is not found in the experimental data. From the strong steady 
state shock indicated by the pressure distribution in figure 22, one should expect from 
past experience a fairly large pressure pulse in the unsteady solution. 

Figures 34 to 39 present for M = 0.8 the jump in unsteady pressure across a NACA 
64A010 airfoil oscillating in plunge at an amplitude of 0.05 chord for reduced 
frequencies of k = 0.05, 0.101, and 0.151. In these examples, the nonlinear time 
dependent solution of EXTRAN2 predicts the shock pulse position more accurately than 
the frequency domain solution. The program EXTRAN2 produces a larger jump in 
pressure in the region between the leading edge and the shock than either the 
frequency domain method or the experiment. This discrepancy is not understood. The 
frequency domain method is in closer agreement with experiment although the 
magnitude for the leading edge is still overpredicted. 

For the NACA 64A010 airfoil oscillating in pitch at M = 0.8 and zero mean angle 
of attack, the theoretical lift and moment coefficients are compared with the 
experimental data of Davis and Malcolm (ref. 14) in figure 40. The real component of 
the lift is in good agreement with experiment, but the imaginary component of the 
theory shows a trend opposite to the experiment in the frequency range 0.05 to 0.2. In 
both theory and experiment, the imaginary component of the lift goes to zero as the 
frequency goes to zero although the theory drops more sharply than the experiment. 
Figure 40 also includes data of Isogai (ref. 15) which generally follows the trends of the 
present theory. 



Figure 41 compares the real and imaginary parts of the moment coefficient about 
the leading edge for the same airfoil oscillating in pitch at M = 0.8 for several values of 
reduced frequency. The real component of the moment from the theory is greater than 
the experimental values but appears to have nearly the same trend. As in the lift case, 
the imaginary component goes to zero more rapidly than the experimental data. The 
results of Isogai (ref. 15) are generally in accord with the present theory. 
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7.0 FLUTTER ANALYSES OF A NACA 64A010 AIRFOIL 


Typical section, two-degree-of-freedom flutter analyses of a NACA 64A010 airfoil 
are presented and discussed in this section. The flutter results were calculated using the 
harmonic aerodynamic procedures discussed in sections 4 and 5. Results from two of the 
three configurations studied are compared with results obtained by Isogai in 
reference 15, using time-integration procedures. Generally, the flutter characteristics 
determined from the two aerodynamic procedures compare favorably. Flutter analyses 
using the time-integration methods predict a deeper flutter bucket than the analyses 
using harmonic procedures. There are no experimental results available for correlation 
with the analytical results of this section. 

The reader is referred to papers by Ashley (ref. 1), Isogai (ref. 15) and Yang et al 
(refs. 31 and 32) for comprehensive reviews of papers and reports concerned with flutter 
analyses of airfoils in transonic flow and only a few comments will be offered here to 
put the current investigations in perspective. First, there have been questions raised 
concerning the flutter results published by Traci et al (ref. 26) and whether or not these 
results include the effects of shock motion. Section 5 of this report discusses this point 
in detail, and it is now hoped that flutter results obtained with harmonic procedures 
will be recognized as true transonic flutter analyses, which include the effects of shock 
motion and the resulting pressure pulse. 

While everyone agrees that transonic flow is a nonlinear phenomenon, there 
remains controversy as to the extent flutter analyses have to be performed in nonlinear 
fashion to obtain useful results. In this respect, the results of Isogai (ref. 15) are of 
special interest for comparison with calculations of this investigation. Isogai uses a 
nonlinear aerodynamic procedure to calculate time-dependent pressure distributions, 
but extracts only the first harmonic of these pressure distributions for use in a linear 
V-g type flutter solution. Generally, Isogai’s results exhibit the flutter boundary 
variations in the transonic regime that were obtained during the present investigation. 
Isogai also has an interesting discussion of the validity of using linear unsteady 
airforces for transonic flutter analyses. 

The other flutter results, which are relatively comparable to the investigation 
discussed here, are presented by Yang et al in ref. 31. First, it should be noted that the 
NACA 64A010 airfoil studied in reference 31 was scaled down to 8.9% maximum 
thickness-to-chord ratio. Second, both aerodynamic theories used, i.e. LTRAN2 based on 
the time-integration method (ref. 4) and STRANS2/UTRANS2 based on the harmonic 
method (ref. 26), were limited to small reduced frequencies, lower Mach numbers and 
small values of A.i = o>M/(l - M2). Although the results presented by Yang, et al., do 
show a degradation in flutter velocity in the transonic regime, these results do not 
manifest the change in flutter modes that will be discussed below. 
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7.1 EQUATIONS OF MOTION 


The flutter problem is formulated in conventional V-g fashion (see, e.g., Chapter 9 
of Bisplinghoff, Ashley and Halfman, ref. 33). Starting with the sketch showing the 
parameters and sign conventions for a binary system (fig. 42), the equations of motion 
are written as 

mh + SqjQ: + Kj^h = - L 


Sojh + Iqi“ + Kq,o! = My 


(7.1) 


where L and My are the total lift (positive up) and pitching moment (about the elastic 
axis, positive nose up) per unit span. The inertial quantities M, Sa and Iqt, are the mass, 
static unbalance, and moment of inertia per unit span about the elastic axis. Equation 
(7.1) may be rewritten as 


.. O L 

h/b+ x^a + covi'^h/b= 

“ mb 


M, 


9 •• 9 9 y 

Xc^b/b + V a + oiorra a = — 7 

mb-^ 


(7.2) 


where Xq. = So;/(mb), ra = Vl^/(mb2), o>h = VKh/m and coq. = Vk^I^. Assuming 
simple harmonic motion, h = hei«"t and a = and writing the lift and moment in 

coefficient form, reduces equation (7.2) to 


- h - x^a + ^ h = - 


pU^c 


-x^h-ij-a + (oi^lu>Ya = 




2mbo)^ 

2 _ pU^c^ 

^ r\i “ — 


o u2 2 

2mb‘^co'^ 


(7.3) 


when p is the free stream density and U is free stream velocity. The lift and moment 
coefficients may be written as 

Chh _ 

Cl= — • h/b + CL£^ - a 




-Mh ^ 


h/b + CMc,-0! 


leading to the eigenvalue equation written in matrix form as 


[M] +■ 


1 


Trpk" 


[A] -x[K] 




(7.4) 
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where 


[M] 


1 



'a 



[K] 

[A] 


("hK 

f 0 

0 


” ^Lh 

— 

2 

■CLq! 

_%h 



and the unknown eigenvalue X is equal to (o>q;/(w)2 (1 + ig) where g is the structural 
damping. The parameters and k are the mass ratio and reduced frequency, 
respectively. 


The aerodynamic coefficients CLh> CLa, CMh and CMa for an arbitrary elastic axis 
position may be obtained from those calculated for an alternate axis a' using the 
following transformations: 


'Clu ' 


“ 1 0 

0 

0~ 





(a' - a)/2 1 

0 

0 


1 CLa 


> = 

-(a - a)/2 0 

1 

0 


r' 

‘-Mh 

k ^Ma > 


_-(a' - a)^/4 -(a'- a)/2 

(a' - a)/2 

1_ 




The airforces for this report are calculated with an elastic axis at the quarter-chord. 


It should be noted that the airforces are calculated for the measured and modified 
NACA 64010 airfoil described in section 6.1 which is slightly nonsymmetric and 10.6% 
thick. 


7.2 A NOMINAL CONFIGURATION 

Flutter results were first obtained for a NACA 64A010 airfoil with the following 
values for the stiffness and mass parameters: 
x^ = 0.25 

r« =0-5 

coh/^o; ~ 0-3, 0.6 
/i = 25 to 300 

Results for a frequency ratio, of 0-3 are presented in figures 43 through 54. 

Figure 43 shows the flutter velocity coefficient and frequency versus mass ratio for 
purely subsonic conditions: an airfoil of vanishing thickness (i.e. a "flat plate”) at 


M = 0.85 and the airfoil configuration at M = 0.75 and 0.80. These curves are 
considered typical for the subsonic condition. The curves are smooth with both the 
flutter velocity coefficient and frequency decreasing with increasing (i.e. decreasing 
air density). Figure 44 shows the amplitude ratio and phase difference of the plunging 
and pitching motions at the elastic axis versus mass ratio for the same subsonic 
conditions. Here, the amplitude ratio increases (becomes more predominantly plunge) 
and the phase angle decreases with increasing yu. Thus, as the aerodynamic damping 
diminishes, the mode of instability tends towards the first natural frequency at the 
system. 

Figure 45 shows the flutter velocity coefficient and frequency versus mass ratio for 
Mach numbers where the flow is transonic (i.e. mixed subsonic and supersonic flow over 
the airfoil) and significant shocks exist on the airfoil. The flutter velocity coefficient has 
the same general behavior as for the subsonic condition, but with much sharper 
gradients at the lower values of /x. Also, the flutter frequency plot indicates a change in 
flutter mode between M = 0.83 and 0.842. At higher values of ix, the mode of instability 
tends to the first natural mode for values of Mach numbers of 0.83 and below, and to 
the second natural mode at M = 0.842 and above. This is reflected in the plots of 
frequency ratio (fig. 45) and of amplitude and phase angle (fig. 46). Note that there is 
another change in flutter mode in the neighborhood of the Mach number at which the 
shock leaves the trailing edge, a phenomenon that will be discussed with later 
examples. 

Mixed subsonic-supersonic flow appears about M = 0.78. The shock first appears 
just ahead of midchord, and moves aft as the Mach number is increased until it finally 
moves off the trailing edge between M = 0.885 and 0.90. The location of the shock on 
the upper surfaces for M = 0.75 to M = 0.90 may be determined from the static pressure 
distributions shown in figure 47. For even higher Mach numbers, a detached shock 
wave appears in front of the leading edge, and this shock again moves aft as Mach 
number is increased until it becomes nearly attached to the leading edge. The shock, 
and its subsequent motion in conjunction with the motion of the airfoil, results in the 
now well recognized pressure pulse in unsteady pressure distributions. The 
corresponding change in the generalized air forces from the flat plate to the airfoil 
solution may be profound as shown in figures 48 to 51. For this case, the imaginary part 
of the Cma term changes sign (becomes positive) between M = 0.82 and M = 0.88. This 
phenomenon helps to explain the switch to a predominantly torsion mode instability 
described in the preceeding paragraph. The gradients in the generalized forces in this 
Mach number range are severe, and more evaluations are required to properly define 
the behavior in the M = 0.85 to M = 0.9 range. Note that the points are connected 
assuming relatively little change in the air forces above M = 0.90. 

Plots of flutter velocity coefficient and frequency versus Mach number for ix= 50 
8u-e shown in figure 52. The flutter boundary for an airfoil of vanishing thickness is 
shown by the dot-dash line. The corresponding boundary for the NACA 64A010 airfoil is 
shown by the solid line. The change in mode shape between M = 0.83 and M = 0.842 
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results, at this mass ratio, in a narrow "chimney” in the boundary, and a high flutter 
speed over a very small range of Mach numbers. The flutter velocity coefficient for the 
airfoil immediately below M = 0.83 is very slightly less than for the flat plate. The 
flutter boundary for the airfoil at Mach numbers above M = 0.842 is well above the flat 
plate boundary. Thus, the airfoil for /jl= 50 exhibits a slight "transonic bucket” in Mach 
number range M = 0.75 to M= 0.83, and undergoes a change in flutter mode between 
M = 0.83 and 0.842 with a resulting rise in the flutter boundary. Flutter solutions were 
also run for M = 0.885 and M = 0.9; however, for the k values studied, there were no 
flutter crossings. 

As discussed in section 6, the harmonic air forces used in the flutter analyses were 
calculated using two different steady velocity potential distributions; the first assumed 
nonviscous flow, the second used a "viscous-wedge” - a moving wedge-nosed ramp to 
approximate the effect of viscosity on the flow (see ref. 20). The flutter boundary 
calculated, including the viscous ramp effects in the steady state potential, are shown as 
a dashed line in figure 52 to provide an indication of the sensitivity of the flutter 
results to the inclusion of viscosity effects. The effect of the viscous ramp is to move the 
shock slightly forward which results in moving the abrupt rise in flutter velocity to a 
slightly higher Mach number. In this particular case, the results at M = 0.83 are very 
sensitive to the viscous effects, and somewhat sensitive at M = 0.87. The calculations at 
M = 0.885 and M = 0.9 were made only without the viscous ramp. 

Similar results are presented in figures 53 and 54 for /x values of 100 and 300. As 
jji is increased, the transonic bucket at Mach numbers below 0.83 increases slightly, and 
the boundary for M = 0.842 and above drops significantly until it falls below the flat 
plate results. Thus, at = 300, the flutter boundary for the airfoil is below that for flat 
plate except in the immediate vicinity of M = 0.83 to M = 0.842. For this mass ratio, a 
flutter crossing was obtained at M = 0.885. No crossings were obtained at M = 0,9 for 
the range of k values studied. 

The next set of figures (figs. 55 to 63) examines the same configuration for a 
frequency ratio of 0.6. Generally, the flutter boundary patterns are the same as for 
coh/<^a “ 0.3. Thus, figures 55 and 56 show the results for purely subsonic cases, and 
both the flutter velocity and frequency behave as discussed for figures 43 and 44. 
Figures 57 and 58 compare an essentially subsonic case (M = 0.83, with viscous ramp) 
with a transonic case (M = 0.842, with viscous ramp). The results at M = 0.842, show 
two flutter modes. The first, which is critical at all but very low mass ratios, is the 
predominantly plunge mode associated with the first natural mode of the system. The 
second, which is critical only at very low mass ratios is the predominantly torsion mode 
associated with the second natural mode of the system. This phenomenon is 
demonstrated for a second flight condition in figures 59 and 60. Here the results for the 
airfoil without viscous effects are compared with results for the flat plate. 

Plots of flutter velocity coefficient and frequency versus Mach number for 
(oYi/cDa = 0.6 are presented in figures 61 through 63. The flutter boundary 
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characteristics are similar to those for a frequency ratio of 0.3. However, the flutter 
boundary for this case is triple-valued for some Mach numbers. The behavior of the V-g 
curves as the flutter boundary starts to rise is sketched in figure 64. As the Mach 
number is increased, the V-g curve bends over, resulting in two crossings of the g = 0 
axis. As the Mach number is increased further, the excursion of V-g curve into the 
positive g half of the V-g plane is reduced until the curve no longer crosses. 

The result is that the portion of boundary in the neighborhood of Me to Md 
represents lightly damped instabilities. As /x increases, the sudden rise in the flutter 
boundary of the basic mode moves to higher Mach numbers, and the depth of the 
associated transonic bucket increases slightly. Also, as fjL increases the flutter boundary 
associated with the systems second natural mode moves to lower values. Again, 
solutions were run at M = 0.885 and M = 0.9. For M = 0.885, a crossing was obtained 
only at /X = 300 for the k values studied. Crossings were obtained for all /x values at 
M = 0.9. The instability mode at M = 0.9 appears to be tending towards the first 
natural mode of the system again and thus represents a change in flutter mode from the 
Mach number range 0.842 to 0.87. Interpretation of the results at M = 0.885 and 
fi = 300 is not clear, and in connecting the points in figure 63 we have tried to match 
modal amplitudes and phase angle differences. 

7.3 TWO EXAMPLES FOR CORRELATION WITH OTHER ANALYSES 

Isogai, in reference 15, computed two flutter examples using air forces from a 
time-integration scheme. These examples were analyzed with a linear flutter analysis 
using the first harmonic components from air forces calculated with a nonlinear 
transonic small perturbation procedure. Results of reanalyzing these two examples 
using the transonic aerodynamic procedures described in sections 4 and 5 are discussed 
in this section. 

The first example, called "Case A’’ by Isogai, was developed to simulate "the 
vibrational characteristics of a typical chordwise section of a swept wing.” The following 
values are used for the mass and stiffness parameters of a binary system: 
a = -2.0, Xqj= 1.8, r^= 1.865, 1.0, ju = 60 

Plots of flutter velocity coefficient and frequency versus Mach number for "Case 
A” are shown in figure 65. Four results are presented: (1) flat plate, (2) the results of 
reference 15 as estimated from the plots therein, and results using the harmonic 
procedures of this report both (3) without viscous effects and (4) with viscous effects. 

The flutter boundary for this example exhibits the same characteristics as that of 
the proceeding example with the abrupt change in flutter mode in the neighborhood of 
M = 0.85. This example has a significantly larger transonic bucket than the previous 
example, with Up/Cbaioj V^l) for the airfoil at M = 0.84 being half the flat plate value. 
The pattern for boundary based on the harmonic procedure looks quite similar to that 
from the time-integration study. However, the dip for the harmonic solution is not quite 
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as great as that from the time-integration solution. The rise in the flutter speed at the 
initial change in flutter mode appears to occur at the same Mach number for the two 
calculations without viscous effects. The inclusion of viscous effects (which results in 
moving the shock forward on the airfoil) moves the velocity rise to a slightly higher 
Mach number. The change in flutter mode as measured in terms of cof/coq. is 
significantly greater for the harmonic solution than for the time-integration solution. 
Also, the harmonic solution again shows three distinct parts to the flutter boundary. 

The first, at Mach numbers for which the flow is purely subsonic, is the classical flutter 
mode which in this case is mostly plunge. The second, at Mach numbers for which there 
are shocks over the chord of the wing, is a predominantly pitch mode. The third, taking 
over at Mach numbers for which the shock approaches the trailing edge, is again 
related to the predominantly plunge mode. 

The level of the flutter velocity coefficient after the change in mode shape is 
significantly higher for the time-integration air forces than for the harmonic air forces. 
It would be expected that the flutter characteristics at M = 0.9 and M = 1.01 would be 
similar since the airforces for these two Mach numbers are nearly the same (see figs. 
12g, 12h and 14 of ref. 15). 

The second example examined by Isogai, "Case B,” uses the following structural 
parameters: 

a = - 0.30, = 0.5, = 0.7, o^h/^o: ^ =60 

which represents a nominal unswept wing. Plots of flutter velocity coefficient and 
frequency versus Mach number are presented in figure 66. Again, results are presented 
for the flat plate and the harmonic procedures, with and without viscous effects, and are 
compared with the results of reference 15. Again, there is an abrupt change in flutter 
mode shape in the neighborhood of M = 0.83. The transonic bucket in the flutter 
boundary just below M = 0.83 is significantly less from the harmonic analyses than 
from the time-integration analysis. Note that for Case A, the two buckets were very 
similar. The flutter velocity coefficient at M = 0.85, which is just above the Mach 
number at which the flutter mode has changed, matches very well between the 
harmonic and time-integration analyses. 
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8.0 IMPROVEMENTS IN NUMERCIAL PROCEDURES 


Two new numerical procedures were investigated. One method, called the two step 
method, has the possibility of considerable improvement in computing time. The 
difference in potential between the airfoil and a flat plate can be computed with good 
accuracy using a solution region which is smaller than that required for calculating a 
full transonic potential. Even though two calculations are required, considerable 
computer time is saved because the more complicated mixed flow problem is computed 
using a much reduced mesh. 

The second numerical procedure is an alternating direction implicit (ADI) 
technique similar to the relaxation method but does not have the limitation to 
frequencies below the critical value related to Mach number and mesh size. This was 
tested for two-dimensional flow but has more practical application to three-dimensional 
flow where the direct solution has memory requirements beyond the capacity of current 
computing machines. 


8.1 THE TWO STEP METHOD 

On comparing the potential distributions for a harmonic flat plate solution with 
that for the full transonic flow over an airfoil at the same Mach number and reduced 
frequency, it was discovered that the difference in potentials was significant only for 
the region close to the airfoil. Hence, if one could obtain a linearized solution for the 
flat plate potential at a greatly reduced computing cost, then the transonic harmonic 
solution for the difference between the flat plate potential and the potential for the 
transonic flow over the airfoil could be computed for a much reduced mesh region near 
the airfoil. 

The complex potential for mixed flow is assumed to be the sum of the flat plate 
potential ijfi and a potential cp2 ;i.e., 

= i/zj +^2 (8.1.1) 

Here satisfies the linearized differential equation for unsteady subsonic flow or 

Ki//i +1^1 -2icoi^i /e + (co^/e)i/'i =0 (8.1.2) 

^xx ^ 

Substituting equation (8.1.1) in the transonic differential equation for <pi and 
using equation (8.1.2.) lead to the same differential equation for <p2 as that for but 
with source terms which depend upon <pQ and The solution of <^^>2 is carried out in the 
same way as for since the source terms merely introduce additional non-zero right 
hand sides into the system of difference equations. The complete derivation of the 
equations is presented in appendix C. 
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Figures 67 to 68 compare the jump in pressure coefficients for the full perturbation 
potential <pi, for the flat plate solution i/^i and for the difference potential (p2- 
Computations were made on the full mesh. Adding the results of the potentials ^f/i and 
<P2 can be seen to yield the values for <pi. 

To obtain a measure of the improvement in using the solution of the flat plate over 
the full mesh to obtain the complete <pi for the reduced mesh, the transonic solution was 
found for the full mesh and for the reduced mesh in a single step. Figures 69 and 70 
indicate that the accuracy of the solution is not adequate with the reduced mesh. 

The solution for <p2 with the source terms calculated from the flat plate solution 
for the complete 68 by 50 mesh was obtained using a reduced mesh of 51 by 22 extracted 
from the complete 68 by 50 mesh. The combined pressure distribution for and <p2 is 
compared with that for <pi over the complete mesh in figures 71 and 72. The accuracy is 
seen to be satisfactory. 

Even though the flat plate solution was obtained from the direct solution of the 
finite difference equation instead of the kernel function method, considerable reduction 
was found in the cost of obtaining the complete transonic solution. The computing times 
are roughly proportional to the square of the number of grid points. Since the flat plate 
solution is antisymmetric, only one-half of the grid points need be used. Thus the flat 
plate solution requires only about 1/4 of the time for the complete <pi solution. The 
reduced grid of 51 by 22 has only 1/3 the number of points and hence requires only 
about 1/9 of the time for the full solution. This is a combined time of only 36% of the 
original full solution. The time for obtaining the flat plate solution could be possibly 
reduced by an order of magnitude by taking advantage of its pure elliptic form. Thus 
the solution cost could be substantially reduced by using the two step method. 

8.2 THE ADI METHOD OF HARMONIC SOLUTION BASED ON TIME 

INTEGRATION 

Because of the size of the matrix for three-dimensional flow, a preferred method of 
solution is some form of relaxation in which the direct solution of the large system of 
difference equations is obtained by a sequence of smaller steps. Conventional relaxation 
methods for both two dimensions and three dimensions have been found to be limited in 
the range of frequencies for which the method converges (see references 9 and 10). 
However, when the equation with the time variable, t, is retained, solutions of harmonic 
motion may be calculated for any frequency range using relaxation type procedures. 
Under the assumption of small perturbations near the speed of sound, the differential 
equation for the perturbation potential is of the form 

[k - (-T+ D* + *yy - (2*xt + *tt) /' = ° (8-2 D 

Rizzetta and Chin (ref. 19) developed an ADI method for solving this equation 
which is an extension of the method of Ballhaus and Goorjian for small frequencies. 



Equation (8.2.1) is linearized by separating the potential into steady state and 
unsteady potentials under the assumption of small oscillations. The steady state 
potential satisfies the classic nonlinear transonic small perturbation differential 
equation with the time derivative term in equation (8.2.1) deleted. The linear 
differential equation for the unsteady perturbation potential \{fi then becomes 

where u = K — (7 + Ds^^Ox- Applying the operator factorization procedure similar to 
that of Rizzetta and Chin (ref. 19) we can solve equation (8.2.2) for each time step by an 
X sweep given by the difference equation 

25x(^ - /(eAt)= [5x(u5x^^) + 5x / 2 + 5yyi//" (8.2.3) 

followed by a y sweep of the form 

( /(eAt^)+ 25^ - <//) /(eAt)= 5yy -<//”)/ 2 (8.2.3) 

where = i/f(nAt) and the tilda denotes an intermediate value of ip between n and n+ 1 . 
Here 8 x on the left-hand side is the backward difference. The differences on the 
right-hand side are central differences at all elliptic points. However, 8 x(u 8 x) for 
hyperbolic points is the backward difference. 

The linearized boundary conditions that the flow velocity on the wing equal the 
velocity of motion on the surface is 

*y=fx + ft 

We introduce harmonic motion of the airfoil and let 
^n^^ngicon At 

then 

= fx 

and the quantity <p^ approaches constant values at each point as t“>oo. 

The complete derivation of this method is included in appendix D. The boundary 
conditions used on the wake are those applied in the direct solution and in the 
relaxation procedure. 

For the first calculations, the same boundary conditions used in the direct 
harmonic solution were applied to the mesh boundaries, but this introduced reflections 
which caused oscillations in the complex potential (p. When the mesh boundary 
conditions used by Kwak (ref. 34) were applied to the difference equations (8.2.3) and 
( 8 . 2 . 4 ), the oscillations were eliminated allowing <p to converge to complex constants at 
every point of the grid. The simpler first order mesh boundary conditions of Engquist 
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and Majda (ref. 22) were also tried. These boundary conditions for pure harmonic motion 
reduce to those used in the direct harmonic solution and are simpler than the boundary 
conditions of Kwak. They also eliminated the boundary reflections. 

A stability analysis of the difference equation has not been made because of the 
complexity of the three time levels introduced by the second time derivative. A Von 
Neumann analysis of the small frequency form (<ptt deleted) is simple and is presented 
in appendix D. This form of the equation is equivalent to the method of Traci, Albano, 
and Farr (refs. 26 and 27) which has the same frequency limitation as our method. The 
difference equation resulting from dropping the first term in parentheses of equation 
(8.2.4) was found to be unconditionally stable. Since the method of Rizzetta and Chin 
(ref. 19) was found to be stable for all frequencies tried, it appears that equations (8.2.3) 
and (8.2.4) are also convergent for all frequencies. 

Solutions were found for M = 0.9 and reduced frequencies of 0.45 and 0.6, 
which are well above the critical frequency of the mesh. The results for the flat plate 
are compared with the method of Rowe et al. (refs. 35 and 36) in figures 73 and 74. The 
comparison with linearized theory is not quite as good as that for the direct solution. 

These results are preliminary. Choosing the size of the time step to speed up 
convergence and studies to improve the accuracy are needed. A solution for mixed flow 
has not been tried but this should cause no difficulty, as solution convergence is 
independent of whether the flow is subsonic or truly transonic. 


35 


I 



9.0 CONCLUSIONS 


This investigation has centered on further development and evaluation of a 
harmonic procedure for analyzing unsteady transonic air forces. It has shown that the 
moving shock associated with an oscillating airfoil is represented in the pressure 
distributions and generalized air forces as calculated with harmonic procedures. 
Pressure distributions of an NACA 64A010 airfoil have been correlated with 
distributions calculated with both harmonic and time-integration procedures. The 
results appear satisfactory as might be expected for an inviscid solution. 

Two-degree-of-freedom flutter analysis of the NACA 64A010 airfoil have been 
performed. It is found that the flutter boundary through the Mach number region for 
which a shock is attached to the airfoil surface may include a flutter '^bucket” (the 
boundary may be well below that for a flat plate) and a change in mode shape resulting 
in a region of low damping and/or a significant rise in flutter velocity. It would be 
helpful to have experimental transonic flutter results to compare with the analytical 
results. 

Theoretical developments from this investigation include development of an ADI 
method which appears to provide a practical relaxation solution approach for the 
three-dimensional problem. Also, preliminary results are presented from a two step 
procedure in which the full velocity potential is found as the sum of the potential for a 
flat plate and an increment in potential due to the thickness distribution. Both of these 
procedures are promising but require further development and assessment. 


36 



I 


APPENDIX A 

DERIVATION OF THE SHOCK RELATIONS FOR A MOVING 
SHOCK REFERRED TO THE STEADY STATE SHOCK LOCATION 


The partial differential equation for transonic small disturbance unsteady flow is 
given in reference 8 as 

[k - (T+ l)^xj *^xx ^yy “ ^xt ' (*^x *^t)t ^ (A-1) 

Here we have dropped the term <pt <pxx to be consistent with Hafez, Riszk, and 
Murman (ref. 28), whose approach we shall follow in the subsequent analysis. (See also 
Williams (ref. 29) and Seebass, et al. (ref, 30) for similar analyses.) The term <p^ <pxx has 
been found to have only a small effect on the solution. To develop the shock conditions 
we express this equation in conservation form. Equation (A-1) becomes 

= ° (A-2) 

Applying the divergence theorem about a thin region containing the shock yields 

■ W"x A + [kv’x -(')'+ lVx^/2] nx + ['Py]ny-[v?t + '^x]*^tA"° 

where (nt,nx,ny) is the normal to the shock and [ ] denotes the jump in the quantity 
across the shock. 

The continuity of tangential velocity is given by 

[ V X n = 0 ( A-4) 

M = 0 (A-5) 

Equation (A-4) when expanded gives 

["x] [%] Pt] 

"x 'V "t 

or 

[<Px]ny- [<^y]nx = 0 


(A-6) 

(A-7) 

(A-8) 
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Let f(x,y,t) = 0 be the shock surface, then since Vf ~ n, we have 


^yl "x ^ fy/ ^ 

/"x = ft Ax = - 

To simplify equation (A-3) we divide by After substituting 

/2= 1/2 + ('"x'^-v) 

=■ <'’x>[l’x] 


where <v?x> = (^x"^ + 9 x“)^ 2 , we obtain for equation (A-3) 

[n] 


A^x] 


+ <K-(7+ l)'Px>+ r i ~ ~r^ = ° 


We now eliminate [<p] by means of equations (A-7) and (A-8). Thus 
-2m ^ 


-2rif /”v\^ /”t\" 

+<K-(7+1)</’A+ — - 

e"x Vx/ Vx/ 


e = 0 


Substituting equations (A-9) and (A- 10) leads to 

2 dx / . /dx\^ 1 /dx\^ 

--.<K-(7.ivx>-y =0 


^ = <Pq + 3. Real I j | 
X =Xo + a Real j XiC^^^fj 


(A-9) 

(A-10) 


(A-11) 


(A-12) 

(A-13) 


where (pQ is the steady state potential, Xq is the steady shock location and a is the 
amplitude of motion. Expanding about the steady state shock location yields for 
equation (A-5) 

M = [^] +aXi[(/^x]ei^t + ... 


= |Vo ^^1®^^^ <^0 + — 

„ X X _J 
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Setting the coefficient of each power of a equal to zero yields for the first two 
terms in the expansion 


[ll>o]=0 

M [*>0 J - 0 

Similarly, expanding equation (A- 11) yields 

2ai6jXje^^Y e + - (7+ 1)<^Q -(7+l)aipj + 

/dxnX^ 

+ <- (7+ l)v50xx “ 


(A-14) 

(A-15) 


‘XX 


/dxn\/dxi\ . . . 

+ 2a \~~] {— — ) + o(a^) = 

\dy/ \dy j ' ' 


Equating the coefficient of each power of a equal to zero leads to the following 
shock relations 

<^K-(7+ l)>Po^^+(dxo/dy)2 =0 

/dxo\ /dxA 

2i.Xi / e - (7 . 1 ^ 2 ( — I = 0 


(A-16) 


In two-dimensional flow the shock is nearly normal and we obtain 


>=0 

(A-17) 


(A-18) 

•Xl-(7+l)<^l^>=0 

(A-19) 


To apply the shock conditions of equation (A- 19) we need to express equation 
(A-19) in terms of <pQ and v?i. Thus eliminating Xi by equation (A-15) we obtain 

^2icc /e - (7+ ^ 

Now 

(7+1)i^q =K-u 

X 
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and 


Equation (A-20) then takes the form 
<(2ico /e + Ux> 

From figure A-1 on the following page we define 
Jl = [uJ = (ui^+ij-uy) 


(A-21) 


(A-22) 


and 


/^2=<^x>=l/2i 


\j - \-l j 


^ig+2j ■ '^ig+l j 


^i -1/2 - ^i„-3/2 \+3l2 - \+\ 12 


“ig j ' ^ig- 1 j j " ^^ig+ 1 j 

^2 = + 


X; - X: 


io-2 


^i+2 


- X: 


^ S “ 

The shock boundary condition in equation (A-21) then becomes 


(^2 + 2iw/^)('^igj-Vlj) 


^1 




"S*' s 
— + 


+1 ■ ""i, 


VI 'V 


> =0 (A-23) 


This equation replaces the finite difference representation of the differential 
equation when shock boundary conditions are employed. 
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APPENDIX B 

DERIVATION OF THE SHOCK POINT OPERATOR 


When a rapidly decelerating flow is supersonic upstream of the point 
ij (ui-1/2 J the flow becomes subsonic downstream of the point 

hi (ui+l/2J ^ 0), a shock wave then lies between the points i and i + l.To satisfy the 
appropriate jump conditions across the shock, the difference operator for the point i,j 
must be conservative. To obtain such an operator we apply the divergence theorem 


^ 7 • F dv = 


F • n dc, 


to the differential operator expressed in conservation form for the control volume 
consisting of lines drawn midway between consecutive columns and rows of mesh points 
as shown in figure B-1 on the following page. Here n is the outward normal to the 
closed surface. We shall consider only two-dimensional flow, but the generalization to 
three-dimensional flow requires only the addition of the central difference operator for 
<pl^^ at the point i,j,k. 


The basic partial differential equation, equation (17) of Ehlers (ref. 8), has the 


form 

^u«P| -2ico</?jye^ +q<P|=0 

where u = K- (y + D^Ox* Hence the vector F is 
F=uv^ 2 -2i6o^Pj/e, 

X y 

and the approximate evaluation of the surface integral in equation (B-1) yields 
8y [“i+l/2j'fxi+y2j -“1-3/2 %3pj-2‘“(‘fitl/2J-*’i-3/2|j)/e] 

* %+i/2 ■ '>1 %-l/2 * '‘iJ ° 

where i+1/2 denotes the value of the quantity at the point midway between xi and 
xi + i; and similarly, for the other half integer subscripts. Dividing the equation by 
)2 X 2y puts the equation in difference form: 

['^i+l/2 j ‘^Xi+1/2 j "'^i-3/2 j '^Xi- 3 / 2 j]/®x " 2(iw/e) j^'Z’i+l/2 j '^i-3/2 j]/®x 


(B-1) 


V 


*%+l/2'%-l/2 




JZ + q-. ip- • = 0 
y ^ij nj 


(B-2) 
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Figure B-1. -Control Volume for Shock Point Operator 
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Now 


[-^i+ l/2j-'^i-3/2j]/«x= ['^i+l/2j-‘^i-l/2j+^i-l/2j-‘^i-3/2j] /®x 

= [(>4+1/2 - >^ 1 - 1 / 2 ) + (>^1-1/2 - H- 3 I 2 ) Vl j] / (^i+1/2 - ^i- 3 / 2 ) 

= [(>^i+l - >^i-l ) -^Xij + (>^i - >^ 1 - 2 ) Vl j] / [""i+l ^ ""i ■ ""i-l ■ ^ 4 - 2 ] 

In reference 8 we substituted the central difference operator for <pxij the upwind 
difference operator for obtained 

[■^i+l/2 i-*’i-3/2 j]/«x = DXl ['=2j(''iJ-^i-lj)-‘<li-l (■^i-lj-*’i-2j)] /(DX, tDXj) 

+ DXj [c 1 i , j - »>ij ) + d ,i (^ij - j)i /(dX, + DXj) (B-3) 

where DXi = xi - xi.2» DX 2 = xi + i - xi_i, and cii, dii, and C2i are given in equations 
(19), (20), and (26) of reference 8. This form was chosen to make equation (B-3) equal to 
a weighted average of upwind and central difference for (^ 1 ^. An alternate form, with 
perhaps better justification, is obtained by replacing <^xij and <^xi_ij with their central 
difference forms. In place of equation (B-3) we get 

[^i+ 1/2 r«’i-3/2 j]/«x = DX, + («>i-lj-*’i- 2 j)]/(DXi +DX 2 ) 

+ DX 2 [cii('Pi+ij-»ij) + d,i(v>ij-Vi_,j)] / (DX, +DX 2 ) (B-4) 

A second modified form is found by substituting 

'^i+1/2 j = (^ij+‘^i+ij) / 2 

‘^i-3/2,j=(V^'^i-2j) /2 


We then obtain 

(V’it l/2j-^l-3/2j ) /«x - (n+lj+*’ij-*’i-lj-‘'i- 2 j)/(DXi +DX 2 ) (E-5) 


From the standpoint of simplicity and logic, the second modified form is to be 
preferred. The first modified form is identical to the second when the x grid in the 
vicinity of the shock is uniform. 
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The first bracketed term in equation (B-2) yields for the second x derivative term 

[“i-H/2 '“i-3/2j%3/2 j] 

= Di [2ci_i uj_ 1/2 j (<^ij - -^i-l j ) - 2di-l 3/2 j j “ 

+ D2 [ 2 CP 1 + 1/2 j («Pi+ij - «Pij) - 2 diUi_ 1/2 j - >^i-l j)] (B-6) 

with Di = DXi/(DXi + DX2) and D2 = DX2/(DXi + DX2). The y derivative becomes 

^yy ' (%+l 12 - %j-l 12) /*y ' 

where aj and bj are given in equation (23) of reference 8. 

Finally the difference equation at the shock point ij is obtained by substituting 
either equation (B-3),(B-4), or (B-5) into equation(B-2). After some simplification , we 
obtain, for equation (B-2), 

■ (^j + + ^2^1 + ^2^2 - ^1^3 - Qij / 2)'^ij + bj ^^+1 

= (D1E3 + D1E4 - D2^2)'^i-lj ■ ^ 1^4 '^i-2j ' ^^2^1 '^i+1 j 

where 

El = ciui+i/2j-iwcii/e 
^2^^i^i-l/2j -it^dii/e 
^3 = %\ “i-l/2j +i^dii_i/e 

E4 = bi_i Ui.3/2j + icodij.j/ e (B-9) 

For the first modified shock point operator, equation (B-8) remains unchanged but 
E3 and E4 are changed to 


% " ^i-l ^i-1/2 j ^ 

E4 = di_i ui.3/2j-iajdii.i/e 


(B- 10 ) 


For the second modified shock point operator, the difference equation replacing 
(B-8) and (B- 9 ) is 


aj <^ij_i - (aj + bj - Qij / 2 El + E2 - E3).^ij + bj <Pij+i 
= (E4 + E3 - E2)¥’i_i j - El <^1+1 j + E4 <^i_2j 


(B-11) 
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where 


El = D2 Ci Ui+i/2 j - iw/[e (dXi + DX2)] 

E2 = D2di Ui.1/2 j + (dXi + DX2)] 

E 3 = DiCi.i Ui.1/2 j - iw /[6 (dXi + DX2)] 

E4 = Didi.i Ui.3/2 j + iw/[e (dXi + DX2)] (B- 12 ) 



APPENDIX C 

DERIVATION OF THE DIFFERENCE EQUATIONS FOR THE 

TWO STEP METHOD 


The two step method utilizes the simpler solution of the flat plate to reduce the 
grid size necessary to obtain a transonic harmonic unsteady solution for the flow over 
an airfoil. The differential equation for transonic harmonic unsteady flow is given in 
reference 8 by 

j J^K-(7+ l)<PoJ j =0 (C-1) 

We assume that <pi is a sum of solutions, and <p2y or 


<Pl='I^l+'P2 (C-2) 

where if/i is the flat plate solution which satisfies the partial differential equation 

Kii/j ~(2icj/e)tpi +(u>^le)^i=0 (C-3) 

XX yy ^ 

and <p2 is the difference between the transonic and the flat plate solutions. Substituting 
equation (C-2) into (C-1) and simplifying the resulting equation by means of equation 
(C-3) yields 

/ uv ?2 ) - (2io) / e)</52 +V ’2 - e - ico(7- 1)<^0 l '^2 

\ X X yy L xxj 

= -io3(7- /(7+ 1) 


At x = xi, we note that 

= ^i +l/2 j'^i-l/2j ^ 2 ^i+1/2 j~^i-l/2 j 
^i+l/2“^i-l/2 ^i+l'^i-1 

For practical notation in coding, we used 

^i+l/2 

^i-1/2 


then 




ico(7 - 1 ) 
7+1 


^x’^l 


= 2(K-Ui+ij)ci( - ^lij ) -2(K-Uij)di(i//lij - 

-2ico(7-l)(ui+ij-Uij)0iij/^7+ l)(xi+i -Xi_i)] 
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The difference operator of the left hand side is derived in reference 8, and the variables 
ci and di are also defined there. Since we divided by 2 in constructing the difference 
operator of the left hand side, the right hand side becomes 

RHS(U) = (k - Uj+i j)ci(.//li+lj - 4>uj) - (K - Uij)di 

where C 3 i = a){y - l)/[(y + 1) (Xi+i - Xi-i). 

Since the boundary conditions on the wing are satisfied by the flat plate solution, 
we have 

d^2 -l<x<l, y = 0. 

The jump condition in <p2 at the trailing edge is treated the same way as the full 
solution and the wake boundary conditions are identical. 

The boundary conditions on the upstream x boundaries of the mesh for <p2 are 
obtained from matching the differential operator d/dx + i<wM/(l-M) of the wing with 
that of the flat plate, i.e., 

- icoMt/?| / (1 - M) = j - icjMi// j (1 - M) 

Since <pi - <p2 we have 

^2 ~ icoM^2 / (I ’ M) = 0 

This is the standard outgoing wave boundary condition which we used for the 
complete solution of <pi. Thus the mesh boundary conditions for <^2 f^e same as for 
the full perturbation potential. 

At h 3 rperbolic points, the flat plate solutions satisfy the elliptic difference equation 
while the complete potential satisfies 

- \^2i (^ij - j) j - '^i-2j)] 

+ 5yy<Pij/2 + [coV 2e + C 3 i(u.+ij - Ujj)] = 0 
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Then substituting <p = + <P2 yields 


■‘^ 2 i-lj) “ ‘^i-l^i-lj('^ 2 i-l j " '^ 2 i- 2 j ) '*’ 1 ij / ^ 

- (We)[c2i('^2ij ”'^2i-lj) " ('^2i-lj ‘ '^2i-2j)] 

+ / 2 e + C 3 J j ■ V> 2 j 

- {<^i-l^ij('^lij - ’^li-lj) -di-lUi-lj('^li-lj - ’^li- 2 j) 

- (ico/e) ^2i ,(’^lij ~ '^li -1 j) ' ‘^li -1 ('^li-lj " '^li- 2 j)] 

+ 5 yy^iij /2 + [coV 2 e + C 3 i(ui+ij - j 

Since satisfies the flat plate equation, then 

+ (i^/e)[cii('/'ii+ij - ’/'lij) + dii ('/'lij - ^li-lj)] 

and the right hand side becomes 

RHS(U) = (Kci-icocii /e) (.//ji+lj - .//lij) 

+ [Kdi-Ci_i Ujj-i(oj/e) (dii-C2i)] {'P uj ~ 'I' U-l j) 

+ (di-lUMj-i^dii_i/e) ('/'jm j - <^ii_2j) 

" ^3i (^i+1 j " ) '^lij 
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APPENDIX D 

DERIVATION OF THE ADI METHOD FOR FREQUENCY DOMAIN 
BASED ON TIME INTEGRATION 


To treat the unsteady flow at transonic speeds we separate the potential <p into 
steady and unsteady potentials 

ip(x,y,t) = i^Q(x,y) + i//(x,y,t) (D-1) 

The differential equation for <p is given in reference 8 by equation (12) as 

[k - (T + 1 )v?x] *^XX ^yy ' (^^xt = ^ (^-2) 

where we have dropped the terms ^t<;^xx and (px!^xt which have (been found to have little 
effect on the solution. The steady state potential (pQ satisfies equation (D-2) without the 
last two time derivative terms. Then assuming small unsteady perturbations, we find 
the linear equation for iff on substituting equation (D-1) becomes 

p - (7 + 1 )v5qJ >^xx - + '^yy " + "^tt) ° (D-3) 


or 


(^^x) / ^yy-(2’^xt + ’^tt)/" = 0 

where 


(D-4) 


u = K-(7+ 1)^0^ (D-5) 

The ADI method of Rizzetta and Chin (ref. 19) for solving equation (D-4) is given by two 
sweeps for each time step. Modifying the difference equations for our particular linear 
partial differential equation yields 
For the x sweep: 

/(^^t)= [Sx(^5^^ ) + 5x(u5^.^")] / 2 + Syyi//" (D-6) 

For the y sweep: 

(<,■>+1 -2*" + *"-') /^At>26,(*"+> - j)/(EAt) 

= Syy(*"'"'-*")/2 (D-7) 

For the boundary conditions on the wing, we use 

<//y = fx + ft (D-8) 
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We assume harmonic motion for the airfoil. Then equation (D-8) becomes 


v!/y = (D.9) 

We now introduce = <pneia>nAt and the boundary conditions become 

= f^ + icx^f (D-10) 

Substituting into equation (D6) and factoring lead to the following 

difference equation for the x sweep: 

2V^eAt)-6^(u5^J-)/2 

= jS] [ Sj (u5j,v>") jl* 6yy*>" + 26y»>“/(e4t)] (D-11) 

For the y sweep we obtain 

(»>"+>- 2/3, »>" + /!, 2 25^ -?)/(eAt) 

-5yy(*^+> -/3,/')/2 (D-12) 

where 

FINITE DIFFERENCE FORMULAS FOR THE X SWEEP 

We shall solve equation (D-11) for *!p. Hence we write 

26,?/(e4t)-8^(u5^?)/2=/3,[26^v"/(eAt)+SyyV" + Jx(u6^,>”)/2] (D-13) 

Following Rizzetta and Chin (ref. 19), we shall use a backward or upwind difference for 
the first derivative with respect to x, and for the remaining space derivatives we shall 
use the form programmed in A344 and described in reference 8. Thus we have for 
elliptic points, with 

C4i = 2/[eAt (xi - Xj.i)] 
the difference equation 

C4i(?ij - j) - CjUj+i j - ) + d^n..(jp.. - j) 

= ^1 [‘=4i('^ij-'^Mj) -V^fj) -2bj(v>5-^^-+l) 

^i’^i+1 j('^i+l j " ) ■ *^rtj(‘^ij “ “^i-l j) J (D-14) 
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For coding the coefficients we write the equations in the form 


SUB 1(1) * ?'i_2j + SUB(I) * + DIAG(I) * J'jj + SUPER(I) * = RHS(I) (D-15) 

At elliptic points this becomes 

SUB 1(1) =0.0 
SUB(I) =-C4i-djUij 
DIAG(I) = C4i + cju-+j j + djUjj 

SUPER(I) = -CiUj+ij 

RHS(I) = (S,[c4i(^!] -v>"lj) + 2aj(v>fj_i -2bj(,.n- 

At hyperbolic points we use backward differences for the second derivative terms in x 
and the difference equations becomes 

‘^4i(‘^ij “'^i-lj) " ‘^i-l'^ij('^ij "‘^i-lj) ‘^i-l^i-1 j ('^i-1 j " '^i-2j) 


= ^1 [^4i(‘^rj +2aj(<pg_i -<pfj)-2bj(^g -^g+i) 

S-l^ij ('^ij ~'^i-lj) ■ ^i-l'^i-lj('^i-lj "'^i-2j)] 

The resulting coefficients in equation (D-15) are 

SUBl(I) = -di,iUi_jj 

SUB(I) - “ % + Cj.jUjj + dj_jUj_j 
DIAG(I)= C4j-Gj_jUjj 
SUPER(l) = 0.0 

RHS(l) = (J,[c 4 i(,.!'j +2aj ) - 2bj (,.!] -v>rj+i) 

+ -^Mj) -<2j)] 


(D-17) 


(D-18) 
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DERIVATION OF MESH BOUNDARY CONDITIONS 

For the boundary conditions on the mesh we use the simple nonreflecting 
boundary conditions derived by Engquist and Majda (ref. 22). For the upstream x 
boundary and our particular form of the equation we have 

-M) = 0 

which becomes 

^X2-Xj) -M " ’^l+l/2j)/f^ -M)At]=0 

in which we apply the boundary conditions at x = (xi + x 2 )/ 2 . Since 

*"+1/2= (*"+l + *")/2, *1+1/2 =(*1+ * 2 ) /2 


we have 

n+1 _ _ *n+l 

where 

Xj -^kl^lj) 

^ki = 0 -^ki)/ (1 + ^ki) 

"kl=M(x2-Xi)/[(l-M)At] 


Similarly the lower outgoing wave type boundary condition is 

i//y - M VK /(l - M^) = 0 


and we get 


(*'■+ 1/2 . *;-+l/2)/ . yj) . M VK - *f 3 / 2 ) /[(l - m2)m] 


At y = (yi + y2)/2 this becomes 

-H+1 n+1 , 

'^il -Ck 2 ‘^i 2 +Xi 


(D-19) 


;D-20) 


(D-21) 


(D-22) 


0 


(D-23) 
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where 


\2 - (1 - \2)/ + c^2) 

C]^2 = M Vk ( y 2 - yi)/[( 1 - m 2) At] 

Xi (v?f2-Ck2'^fl) 

For the downstream x boundary, the boundary conditions at 
^ ~ ^^imax *imax‘ 1^^2 are 

+ +M) = 0 

which becomes 

, n+1 _ n+1 

^maxJ ^max'^j 

where 

'k3 = (>-‘^3)/0+‘>k3) 

('"Lx-ij ■'^3 


(D-24) 


(D-25) 


(D-26) 


(D-27) 

(D-28) 


The upper y mesh outgoing plane wave boundary conditions are 

'Py + M\K\p^ l(l -M^) = 0 (D-29) 

In the same manner the difference form of the boundary conditions is seen to be 
,.n+] _ _ n+1 

^jmax ^^4 ^2 (D-30) 

where 

^k4 ~ (* ■ <^k4)/(^ <^k4) (D-31) 

and 


X2=^\ 

^ ^ V ^max 




-l-'^k4^ij^3^ 


(D-32) 
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APPLICATION OF THE MESH BOUNDARY CONDITIONS TO 
THE DIFFERENCE EQUATIONS FOR THE X SWEEP 

For i = 2, we apply the upstream x boundary condition, namely, equation (D-21), 
which becomes 

^lj = ^kl^2j+Xo (D-33) 

Substitution into equation (D-15) yields, with the aid of equation (D-16) 

DIAG(2) = DIAG(2) + SUB(2) (D-34) 

RHS(2) = RHS(2) - Xq SUB(2) 

SUB(2) = 0.0 (D-35) 

where 

X0 = ^l('^2j -^kl'^lj) 

Boundary conditions on the upper and lower boundaries have no effect on the x 
sweep. 

For i = imax — 1, we apply the downstream boundary conditions 

i ~ ^k3 -1 i ^3 (D-36) 

^maxJ ^max ^ 

Then we modify DIAG, SUPER, and RHS by 

DIAG(I) = DIAG(I) + cj^3 SUPER(I) 

RHS(I) = RHS(l) - X3 SUPER(l) 

SUPER(I) = 0.0 (D-37) 

for I = IMAX - 1. Here, 

X3 = ^1 (’^imax-1 j ' ""kS %axj) 

BOUNDARY CONDITIONS ON THE WING 

We apply the boundary conditions on the airfoil at y = 0 which lies halfway 
between the j = jm and the j = jm + 1 y grid points. We apply the boundary conditions 

<fiy = (D-38) 
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where FiL is given by equation (D-10) for the lower surface. In difference form this 
becomes 


h = F.L 


Hence 


where 


and 


-b: +l) = b: hP/ 

bj = l/(yj+, -yj.i) (yj+i-yj) 


(D-39) 


h = Vj +1 

■'m -'m 

Then the right hand side term is modified by 
RHS(I) = RHS(I) + 2f 1 (h FjL + 
Similarly for j = jm + 1 we have 

+1 Wi +l')^"^j +l^^i^ 

Jm ^ V ^Jm ^Jm V Jm ^ ^ 

and the right hand side is modified by 

RHS(I) = RHS(I) - 2/3iaj^+i (h 


(D-40) 


(D-41) 


(D-42) 


We apply the same harmonic boundary conditions in the wake as in the A344 
program. Thus, the jump in potential across the wake is given by 

icoCx-x+.'l 
Ai/5(x) = ^ 

We add A(p to the (pyy difference to make (py continuous across the wake. From equation 
(D-16) we have for j = jm 


^1 ['=4i 

(’’L ■ 

n > 


- ^‘■in. ( 

/ n 

n n 

^ n\ 


-^ij 

„+l+^uj 


/ n 

n 

\1 



)J 


(D-43) 


Defining A<pin = tpijn - (pjn^ we obtain 

RHS(I) = RHS(I) - 2j3 1 b: A</> " 
^ Jm 
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Similarly, for j = jm + 1- 

RHS(I) = RHS(I) + 2|3iaj^+i (D-44) 

The boundary conditions on the wing for h 3 rperbolic points are the same as for 
elliptic points, and all mesh boundary conditions are subsonic. 

DIFFERENCE EQUATIONS FOR THE Y SWEEP 

We solve the difference equation (D-12) for the n+1 time step. Thus we rewrite 
equation (D-12) in the form 

6yy«."-" V 2 - /(eAt2)- /(.At) 

= |3, [Syy*;"/2+ (|S,»’‘'‘'-2*>")/(eAt2| -26j?/(eAt) 

Writing out the differences, we get 

(D-45) 

where C 4 j = 2/[eAt(xi - xf-i)], Ei = l/(eAt2) and y8i = 

In the same manner as in the x sweep we define the coefficients 

SUB(J) * + DIAG(J) >1= + SUPER(J) * ipJ+j = RHS(J) (D-46) 

Then 

SUB(J) =aj 

DIAG(J) = - aj - bj - Ej - C4i 
SUPER(J) = bj 

RHS(J) C4i (?ij +?, [aj (,>!]., 

-bj (*>rj-f"j+i) -2Ei*.|] +^,Ei,pS-‘] -C4iv!itj (D-47) 
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Note that the right-hand side contains a term of the n +1 step. Since we are 
sweeping in the direction of increasing i, this term is known from the previous column 
solution. This sweep replaces the values of the values of Hence the old 

values of must be saved before written over by the n +1 value of prior to 
going to the ith column calculations. 

APPLICATION OF THE MESH BOUNDARY CONDITIONS FOR THE Y SWEEP 

For boundary conditions on upstream boundary ,i + 2, we have from equation 
(D-23) 


n +1 ^n+ 1 _ 

='^2j ‘^kl+Xj 


(D-48) 


where 


Substituting into equations (D-46) and (D-47) modifies the coefficients according to 
DIAG(J) = DIAG(J) + cj^i C 42 

RHS(J) = RHS(J) - C42Xj + 

The last term is to correct for the last term in RHS(J) in equation (D47). 

For the revised coefficients on the lower mesh boundary, j = 2, we substitute the 
boundary conditions 

n +1 ,n+l_ . 

^il ='^i2 ^k2+Xi 

into equation (D-45) and modify the coefficients in equation (D-46) according to 
DIAG(2) = DIAG(2) + Cj^2SUB(2) 

RHS(2) = RHS(2) - SUB(2) Xi D-49) 


where 


SUB(2) = 0.0 


For the revised coefficients on the upper mesh boundary, j = jmax ~ 1. we 
substitute the boundary conditions 

n +1 _ n +1 _ 

*^ii *^ii "I ^k4 Xo 
^Jmax ^Jmax ^ ^ 
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into equation (D-45) and then modify the coefficients in equation (D-47) according to 
DIAG(JMAX - 1) = DIAG(JMAX - 1) + SUPER(JMAX - 1) Cj^4 
RHS(JMAX - 1) = RHS(JMAX - 1) - SUPER(JMAX - 1) X2 
SUPER( JM AX - 1 ) = 0 .0 (D-50) 


where 


The difference equation is not affected by the downstream boundary conditions. 
However, since all mesh boundary values of <p must be computed after each sweep we 
need to use equations (D-25) and (D-27) to find ¥^imaxj- 


For the revised coefficients for the airfoil boundary conditions, we apply the 
boundary conditions of equations (D-39) to equation (D-43). This yields for the second 
term of equations (D-45) at j = jm 




_n+l 



hFjl- 
m ^ 


(D-51) 


This term then contributes to the RHS term with a change of sign and eliminates 
terms from DIAG and SUPER. A term similar to equation (D51) is obtained from the 
right-hand side of equation (D-45). Thus the coefficients in equation (D-47) are modified 
by 

DIAG(JM) = DIAG(JM) + b: 

Jm 

SUPER(JM) = 0.0 


RHS(JM) = RHS(JM) + ^h Fj^ + ) - h J (D-52) 

We now consider the terms for the upper airfoil boundary, j = jm + 1* From 
equation (D-41), we have 


^ +1 




hFj^ 


(D-53) 
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This again contributes an additional term to the right-hand side. Applying 
equation (D-53) the right hand side of equation (D-47) yields 

DIAG(JM + 1 ) = DIAG(JM + 1 ) + a: +i 

■'m ‘■ 

SUB(JM + 1) = 0.0 

RHS(JM + 1) = RHS(JM + 1) - ^h -t- " h J (D-54) 

For j = jm and j = jm + 1, we must also revise the coefficients for the wake boundary 
condition. For j = jm, we obtain for the second term of equation (D-45) 






n+1 






-b| 

Jm ^ 


n+1 


(D-55) 


This adds to the RHS term. Also applying equation (D55) to the similar 

terms of the RHS in equation (D-47), we obtain for the modified right-hand side 

RHS(JM) = RHS(JM) + - /^l (D-56) 

Similarly for j = jm + 1, we get 

RHS(JM+ 1) = RHS(JM+ D-aj^-hl (D-57) 


STABILITY ANALYSIS OF THE METHOD FOR SMALL FREQUENCY 

Because of the three steps in time required by the <ptt term, a von Neumann 
stability analysis of the difference equations (D-11) and (D-12) is difficult and has not 
been made. Dropping the <ptt term is equivalent to assuming small frequency. The 
corresponding frequency domain method was derived by Traci, Albano, and Farr 
(refs. 26, 27) and was solved using the relaxation technique. Their difference equations 
have the same frequency limitation as our equations. Since a stability analysis of the 
equation neglecting the term <^tt is simple we shall analyze the difference equations 
(D-11) and (D-12) with the first term of equation (D-12) deleted, or 

2S,v/(eAt)-6^(uS^?)/2=(Si [2sX/(eAt)+ 6^ (uS^v”)/ 2 + 6yy,>"] (D-58) 
2S^/'^'/(eAt)- Syy'p"'^ V 2 = (3, "/ 2 (D-59) 
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The error in for equally spaced points in the grid may be represented by an 
expression of the form 

E|Jj = 53 exp'[i(k£7rAx/a+jmA/b)] 
fi,m 

and the stability of each component will be examined separately. Now, we have for the 
flat plate, 

25x</’”/ eAt ~ O!^ 1^1 - exp(-i7r£ Ax /a)j (D-60) 

^x j - + <pf_i j) /ax2 ~ - 2oi2<P^ (D-61) 

where 

aj =2/^eAtAx), .a 2 = K ^1 - cos — j Ax^ 

and an exponential factor common to all terms is not written. Similarly 

5yyV?” - loL^'P^ (D-62) 

where «3 = (1 — cos m7rAy/b)/Ay2. With equations (D-60), (D-61), and (D62) substituted 
into equations (D-58) and (D-59), we obtain, with 

^[“1 "^“ 2 ] ■“2-2“3] 

[oij ( 1 - e +“ 2 ] “ “1 (1 " + 

where 6\ = tt Ax/a. 

Eliminating Ip from the second equation by means of the first yields the following 
relation for the amplification factor, after some algebraic manipulations, 

“1 («i(i -ag 

“1 (1 - +"2) (“1 C " ® 

The magnitude of jSi is unity and the magnitude of the two factors in curly 
brackets is seen to be less than unity since they are of the form 

[(a-b)-ic] /[(a + b)-icl 

for a,b,c positive and real. Thus the time dependent method equivalent to the relaxation 
method of Traci, Albano, and Farr (ref. 26) converges for all frequencies. Since the time 
dependent method of Rizzetta and Chin has been found to be stable for a wide range of 
frequencies it appears that the ADI method for harmonic motion described by equation 
(D-11) and (D-12) is also convergent for all frequencies. 


n+1 
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Figure 1.— Instantaneous Unsteady Pressure Distribution on the Upper Surface of a 
Pitching NACA 64A010 Airfoil Computed From the Harmonic Solution 
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near shock 
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Figure 3.— Influence of Grid Refinement on the Rea! Part of the Pressure Pulse, NACA 

64A010 Airfoil Oscillating in Pitch About the Quarter-Chord Point, M=0.8, k-0.3 
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□ Ax = 0.02 at shock 
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Figure 6.— Influence of Grid Refinement on the Jump in imaginary Part of Unsteady Potential 
Across a NACA 64A010 Airfoil Oscillating in Pitch About the Quarter-Chord Point 
M=0.8,k=0.3 
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Figure 7. -Three-Dimensional Surface Representation 
Potential, Coarse Grid. M=0.8,k=0.3, a =1 


of Reai Part of Unsteady 
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Figure 8.-Three-Dimensional Surface Representation of imaginary Part 
of Unsteady Potential, Finer Grid, M=0.8, k-0.3, -1 
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Figure 9.— Three-Dimensional Surface Representation of imaginary Part of 
Unsteady Potential, Finest Grid, M=0.8, k=0.3, 
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Figure 1 0.— Three-Dimensional Surface Representation of Real Part of Unsteady 
Potential for Solution With Coarse Grid and First Form of the Shock 
Point Operator, M=0.8,k=0.3,a^=1° 
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Figure 1 3.— Three-Dimensional Surface Representation of Real Part of Unsteady 
Potential for Solution With Coarse Grid and Third Form of the Shock 
Point Operator, M=0.8, k=0.3, 



77 




1.' 


5.i 




Figure 1 4. — Three-Dimensional Surface Representation of imaginary Part of Unsteady 
Potential for Solution With Coarse Grid and Third Form of the Shock Point 
Operator, M=0.8, k=0.3, a^=1° 



I 



79 



CD Shock boundary conditions □ Shock boundary conditions 

^ Shock point operator ^ Shock point operator 



Figure 16.— Comparison of the Jump in Potential Across the Airfoil for Fitted Shock 
From Shock Boundary Conditions With Results From the Shock Point 
Operator, M=0.8, k=0.3, oc^=1° 
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Figure 18.— Effect of Shock Point Operator on Real Part of Unsteady Potential With 
Weak Steady Shock, M=0.8, k=0.3, ol^= 0°. Fine Grid 
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2 Xq = Shock trajectory 



Assumed shock displacement Pressure variation in point A(P^) 

Figure 19.— Contribution of the Periodicaiiy Moving Shock Wave to the Pressure 
Signai at a Fixed Observation Point. (From Tijdeman, Reference 24} 
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Figure 20-Comparison of Steady State Pressure Distribution on a NACA 64A010 

Airfoil as Computed by EXTRAN2 With NASA-Ames Experiments, M=0.75, 




® NASA-Ames experiments (SI 30) 

TSFOl L, Spreiter scaling 

TSFOl L, Krupp scaling 

EXTRAN2, Spreiter scaling with Riegel's rule. 


Figure 21,— Comparison of Steady State Pressure Distribution on a NACA 64A010 

Airfoil as Computed by Three Methods With NASA Experiments, M=0,8, 
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Figure 22.— Comparison of Steady State Pres 
Airfoil as Computer by EXTRAl 
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Figure 23.— Pressure Coefficient Distributions on NACA 64A010 Airfoil Pitching 1 
About Quarter-Chord Point, M-0.75, k=0.20, 







Figure 25 —Pressure Coefficient Distributions on NACA 64A010 Airfoil Pitching 1° 
About Quarter-Chord Point, M=0.8, k=0. 101, a =0 , Imaginary Part 
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Figure 26, —Pressure Coefficient Distributions on NACA 64A010 Airfoil Pitching 1 
About Quarter-Chord Point, M~0.8, k=0.202, ot^=0 , Rea! Part 
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Figure 27.--Pressure Coefficient Distributions on NACA 64A010 Airfoil Pitching 1 
About Quarter-Chord Point, M=0.8, k=0.202, a^=0°. Imaginary Part 
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Figure 28.— Pressure Coefficient Distributions on NACA 64A010 Airfoil Pitching 1 
About Quarter-Chord Point, M=0.8, k=0.247, =0° , Real Part 





© NASA-Ames experiments (Dl 56) 

Present theory 

Time dependent (EXTRAN2) 


Figure 29.— Pressure Coefficient Distributions on NASA 64A010 Airfoii Pitching 1 
About Quarter-Chord Point, M=0.8, k=0.247, = 0°, imaginary Part 





Figure 30.— Pressure Coefficient Distributions on NACA 64A010 Airfoil Pitching t 
About Quarter-Chord Point, M=0.8, k=0.303, =0°, Rea! Part 
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Figure 31 .—Pressure Coefficient Distributions on NACA 64A0W Airfoil Pitching 1 
About Quarter-Chord Point, M=0.8, k=0.303, a =0°, Imaginary Part 





Figure 32.— Pressure Coefficient Distributions on NACA 64A010 Airfoii Pitching 1 
About Quarter-Chord Point, M=0.842, k=0.202, =0° , Reai Part 
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Figure 33,-Pressure Coefficient Distributions on NACA 64A010 Airfoil Pitching 1° 
About Quarter-Chord Point, M=0.842, k=0.202, a =0°, imaginary Part 
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Figure 34— Pressure Coefficient Distributions on NACA 64A010 Airfoil Oscillating 0.05 
Semichord in Plunge, M=0.8, k=0.05, a™ = 0°, Real Part 
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Figure SS.-^Pressure Coefficient Distributions on NACA 64A010 Airfoil Oscillating 0.05 
Semichord in Plunge, M=0.8, k=0.05, a = 0 ° Imaginary Part 
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Figure 36,— Pressure Coefficient Distributions on NACA 64A010 Airfoil Oscillating 0.05 
Semichord in Plunge, M=0.8, k=0. 101, = 0 ^, Real Part 
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Figure 37.— Pressure Coefficient Distribution on NACA 64A010 Airfoii Osciiiating 0.05 
Semichord in Plunge, M=0.8, k=0. 101, a^ = 0°, Imaginary Part 
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Figure 38.— Pressure Coefficient Distributions on NACA 64A010 Airfoii Oscillating 0.05 
Semichord in Plunge, M=0.8, k=0. 151, a^ = 0°, Peal Part 
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Figure 39.— Pressure Coefficient Distributions on NACA 64A010 Airfoil Oscillating 0.05 
Semichord in Plunge, M=0.8, k=0. 151, = 0 °, Imaginary Part 






Figure 41 .'^Comparison With Experiment of the Calculated Lift Coefficient for a NACA 
64A010 Airfoil Oscillating 1° in Pitch About the Quarter-Chord Point, 

M=0.8, afjj=0° 
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Figure 42,-Definition of Parameters for an Airfoil Oscillating in Pitch and Plunge 


106 







M = 0.85 



50 100 150 200 250 





Figure 44.-F!utter Modes (Amplitude Ratio and Phase Difference) Versus Mass Ratio 
Fiat Plate, coa / oj ™ = 0.3 















O O 


k = 0.150 



Figure 49.— Moment Coefficient Versus Mach Number fora NACA 64A010 
Airfoil Oscillating in Plunge 
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Figure 52.— Flutter Velocity and Frequency Ratio Versus Mach Number for a 
NACA 64A010, = 0.3, g = 50 
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0 With viscous ramp 

© Flat plate 



Figure 53.— Flutter Velocity and Frequency Versus Mach Number fora NACA 64A010 
Airfoil, co/j/cOq, =0.3, ii= 100 
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-igure 54.— Flutter Velocity and Frequency Versus Mach Number for a NACA 
AifoU, co/j/cOq, = 0.3, g = 300 


0.95 


64A010 






Figure 56. —Flutter Modes (j 
Ratio, Subsonic 
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Figure 57.— Flutter Velocity and Frequency Versus Mass Ratio for a NACA 64A010 Airfoil 



(cOf^/ojQ,) = 0.6 

With viscous ramp 




Figure 58.— Flutter Modes (Amplitude Ratio and Phase Difference) Versus Mass 
Ratio for a NACA 64A010 Airfoil 
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Figure 59.— Flutter Velocity and Frequency Versus Mass Ratio for a NACA 64A010 Airfoil 
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Figure 60.— Flutter Modes (Amplitude Ratio and Phase Difference) Versus Mass Ratio 
fora NACA 64A010 Airfoil 
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Figure 62.— Flutter Velocity and Frequency Versus Mach Number for a NACA 64 A0 10 
Airfoil, = 0.6, g= 100 
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Figure 63.— Flutter Velocity and Frequency Versus Mach Number for a NACA 64A010 
Airfoil, =0,6, g = 300 
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Case B, M = 60 
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Figure 73.— Comparison of Pressure Coefficient Distribution From an AD! Solution 
With Linearized Theory for the Flow Over a Fiat Plate Oscillating in 
Pitch About the Leading Edge, M = 0.9, k = 0.45 
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Linear theory (Rowe et al., Ref. 35) 
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Figure 74.— Comparison of Pressure Coefficient Distribution From an AD! Solution 
With Linearized Theory for the Flow Over a Fiat Plate Oscillating in 
Pitch About the Leading Edge, M = 0.9, k = 0.6 
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